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1 Introduction 



When probed at very short wavelengths, QCD is essentially a theory of free 'partons' — quarks and 
gluons — which only scatter off one another through relatively small quantum corrections, that can be 
systematically calculated. At longer wavelengths, of order the size of the proton ~ Ifm = 10~^^m, 
however, we see strongly bound towers of hadron resonances emerge, with string-like potentials build- 
ing up if we try to separate their partonic constituents. Due to our inability to solve strongly coupled 
field theories, QCD is therefore still only partially solved. Nonetheless, all its features, across all 
distance scales, are believed to be encoded in a single one-line formula of alluring simplicity; the 
Lagrangian of QCD. 

The consequence for collider physics is that some parts of QCD can be calculated in terms of 
the fundamental parameters of the Lagrangian, whereas others must be expressed through models or 
functions whose effective parameters are not a priori calculable but which can be constrained by fits to 
data. However, even in the absence of a perturbative expansion, there are still several strong theorems 
which hold, and which can be used to give relations between seemingly different processes. (This 
is, e.g., the reason it makes sense to constrain parton distribution functions in ep collisions and then 
re-use the same ones for pp collisions.) Thus, in the chapters dealing with phenomenological models 
we shall emphasize that the loss of a factorized perturbative expansion is not equivalent to a total loss 
of predictivity. 

An alternative approach would be to give up on calculating QCD altogether and use leptons 
instead. Formally, this amounts to summing inclusively over strong-interaction phenomena, when 
such are present. While such a strategy might succeed in replacing what we do know about QCD 
by "unity", however, even the most adamant "chromophobe" must acknowledge a few basic facts of 
collider physics for the next decade(s): 1) At the Tevatron and LHC, the initial states are unavoidably 
hadrons, and hence, at the very least, well- understood and precise parton distribution functions (PDFs) 
will be required; 2) high precision will mandate calculations to higher orders in perturbation theory, 
which in turn will involve more QCD; 3) the requirement of lepton isolation makes the very definition 
of a lepton depend implicitly on QCD, and 4) the rate of jets that are misreconstructed as leptons in the 
experiment depends explicitly on it. Finally, 5) though many new-physics signals do give observable 
signals in the lepton sector, this is far from guaranteed. It would therefore be unwise not to attempt 
to solve QCD to the best of our ability, the better to prepare ourselves for both the largest possible 
discovery reach and the highest attainable subsequent precision. 

In the following, we shall focus squarely on QCD for mainstream collider physics. This in- 
cludes factorization, hard processes, infrared safety, parton showers and matching, event generators, 
hadronization, and the so-called underlying event. While not covering everything, hopefully these 
topics can also serve at least as stepping stones to more specialized issues that have been left out, such 
as heavy flavours or forward physics, or to topics more tangential to other fields, such as lattice QCD 
or heavy-ion physics. 

1.1 A First Hint of Colour 

Looking for new physics, as we do now at the LHC, it is instructive to consider the story of the 
discovery of colour. The first hint was arguably the A++ baryon, found in 1951 [ ]. The title and 
part of the abstract from this historical paper are reproduced in figure 1 . In the context of the quark 
model — which first had to be developed, successively joining together the notions of spin, isospin, 
strangeness, and the eightfold way — the flavour and spin content of the baryon is: 




(1) 
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Meson- Nucleon Scattering and Nucleon Isobars* 

Keith A. Brueckner 
Deparimeni of Physics, Indiana University, Bloomington, Indiana 
(Received December 17, 1951) 



"f...J It is concluded that the apparently anomalous features of the scattering can be interpreted to 
be an indication of a resonant meson-nucleon interaction corresponding to a nucleon isobar with 
spin |, isotopic spin |, and with an excitation energy of 277 MeV." 

Fig. 1: The title and part of the abstract of the 1951 paper [.] (published in 1952) in which the 
discovery of the baryon was announced. 



clearly a highly symmetric configuration. However, since the A"*""*" is a fermion, it must have an over- 
all antisymmetric wave function. In 1965, fourteen years after its discovery, this was finally under- 
stood by the introduction of colour as a new quantum number associated with the group SU(3) [ , ]. 
The A"*""*" wave function can now be made antisymmetric by arranging its three quarks antisymmetri- 
cally in this new degree of freedom, 

\A++) = e'^''\ui^Uj^Uk^) , (2) 

hence solving the mystery. 

More direct experimental tests of the number of colours were provided first by measurements 
of the decay width of vr" — 77 decays, which is proportional to N^,, and later by the famous "R" 
ratio in e+e" collisions. Below, in Section 1.2 we shall see how to calculate such colour factors. 

1.2 The Lagrangian of QCD 

Quantum Chromodynamics is based on the gauge group SU(3), the Special Unitary group in 3 (com- 
plex) dimensions. In the context of QCD, we represent this group as a set of unitary 3x3 matrices 
with determinant one. This is called the adjoint representation and can be used to represent gluons 
in colour space. Since there are 9 linearly independent unitary complex matrices, one of which has 
determinant —1, there are a total of 8 independent directions in the adjoint colour space, i.e., the glu- 
ons are octets. In QCD, these matrices can operate both on each other (gluon self-interactions) and 
on a set of complex 3-vectors (the fundamental representation), the latter of which represent quarks in 
colour space. The fundamental representation has one linearly independent basis vector per degree of 
SU(3), and hence the quarks are triplets. 

The Lagrangian of QCD is 

C = i,\{ir){D^)i,% - mg^i^iiPq, - , (3) 

where denotes a quark field with colour index i, ipq = [ipqR, ihi/O, ipqBy , 7^ is a Dirac matrix that 
expresses the vector nature of the strong interaction, with ^ being a Lorentz vector index, ruq allows 
for the possibility of non-zero quark masses (induced by the standard Higgs mechanism or similar), 
F^^ is the gluon field strength tensor for a gluon with colour index a (in the adjoint representation, 
i.e., a e [1, . . . , 8]), and is the covariant derivative in QCD, 

{D^,)ij = 5ijd^ - igst^Al , (4) 
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Fig. 2: Illustration of a gg(7 vertex in QCD, before summing/averaging over colours: a gluon in a state 
represented by interacts with quarks in the states ipgR and i/jqc- 



with gs the strong coupling (related to Os by = Anas; we return to the strong coupling in more 
detail below), the gluon field with (adjoint-representation) colour index a, and proportional to 
the hermitean and traceless Gell-Mann matrices of SU(3), 
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These generators are just the SU(3) analogs of the Pauli matrices in SU(2). By convention, the 
constant of proportionality is normally taken to be' 

tij = l^j ■ (6) 

This choice in turn determines the normalization of the coupling gs, via equation (4), and fixes the 
values of the SU(3) Casimirs and structure constants, to which we return below. 

An example of the colour flow for a quark-gluon interaction in colour space is given in fig- 
ure 2. Typically, however, we do not measure colour in the final state — instead we average over 
all possible incoming colours and sum over all possible outgoing ones, wherefore QCD scattering 
amplitudes (squared) in practice always contain sums over quark fields contracted with Gell-Mann 
matrices. These contractions in turn produce traces which yield the colour factors that are associated 
to each QCD process, and which basically count the number of "paths through colour space" that the 
process at hand can take, modulo that the convention choice represented by equation (6) introduces a 
"spurious" factor of 2 for each power of the coupling Og, as we shall see^. 

A very simple example of a colour factor is given by the decay process Z — )• qq. This vertex 
contains a simple 6ij in colour space; the outgoing quark and antiquark must have identical (anti-)col- 



' Another choice that is occasionally (though rarely) seen in the literature is t — A/\/2- This gives a more intuitive 
colour counting, but since it also implies a different normalization for the coupling and since most text material uses the 
convention defined by equation (6), we shall stick to that choice for the remainder of these lectures. 

^Again, although one could in principle absorb that factor into a redefinition of the coupling, effectively redefining the 
normalization of "unit colour charge", the standard definition of as is now so entrenched that alternative choices would be 
counter-productive, at least in the context of a supposedly pedagogical review. 
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ours. Squaring the corresponding matrix element and summing over final-state colours yields a colour 
factor of 

Z^qq : ^ \M\^ ^ 5,j6*i = Tt{6} = Nc = 3 , (7) 

colours 

since i and j are quark (i.e., 3-dimensional fundamental-representation) indices. 

A next-to-simplest example is given by the Drell-Yan process, qq — )• 7*/^, i.e., just a crossing 
of the previous one. By crossing symmetry, the squared matrix element, including the colour factor, is 
exactly the same as before, but since the quarks are here incoming, we must average rather than sum 
over their colours, leading to 

^ ^ ■ In 1^1' ^ l^^^^i^ = = 1, (8) 

colours 

where the colour factor now expresses a suppression which can be interpreted as due to the fact that 
only quarks of matching colours are able to collide and produce a Z boson, effectively reducing the 
incoming quark-antiquark flux by a factor l/Nc- 

To illustrate what happens when we insert (and sum over) quark-gluon vertices, such as the 
one depicted in figure 2, we take the process Z — )• 3jets. The colour factor for this process can be 
computed as follows, with the accompanying illustration showing a corresponding diagram (squared) 
with explicit colour-space indices on each vertex: 

Z ^ qgq : 

colours 

= Iti{6} = 4, 

where the last Tr{5} = 8, since the trace runs over indices in the 8-dimensional adjoint representation. 

The tedious task of taking traces over SU(3) matrices can be greatly alleviated by use of the 
relations given in Table 1. In the standard normalization convention for the SU(3) generators, equa- 
tion (6), the Casimirs of SU(3) appearing in Table 1 are^^ 

Tr = ^ Cf = ^ Ca = Nc = ?>. (10) 

In addition, the gluon self-coupling on the third line in Table 1 involves factors of /"*'^. These are 
called the structure constants of QCD and they enter due to the non-Abelian term in the gluon field 
strength tensor appearing in equation (3), 

= d,Al - d,Al + gsr^'A^Al . (1 1) 

^Se.e, e.g., [5, Appendix A.3] for how to obtain the Casimirs in other normalization conventions. 
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Trace Relation 



Indices 



Occurs in Diagram Squared 



ah 



ab 



a,6e [1,...,8] 



ae [1,...,8] 
i,j,k G [1, ... ,3] 



a, b,c,d ^ [1, 



tf^tl^ = Tr (^SJkS^e - t^S^jSm) i, j, /c, ^ G [1, . . . , 3] 




(Fierz) 



Table 1: Trace relations for t matrices. These relations are convention-independent as they stand. 
Relations for a specific normalization convention for the t matrices are obtained by inserting the 
specific values of Tr, Cr, and Ca pertaining to that convention choice, as discussed in the text. More 
relations can be found in [4, Section 1.2] and in [5, Appendix A.3]. 



The structure constants of SU(3) are listed in 
the table to the right. Expanding the F^yF^^'^ 
term of the Lagrangian using equation (11), we 
see that there is a 3-gluon and a 4-gluon vertex 
that involve /"^'^, the latter of which has two 
powers of / and two powers of the coupling. 
Finally, the last line of Table 1 is not really 
a trace relation but instead a useful so-called 
Fierz transformation. It is often used, for in- 
stance, in shower Monte Carlo applications, 
to assist in mapping between colour flows in 
Nc = 3, in which cross sections and split- 
ting probabilities are calculated, and those in 
Nc — ?• oo, used to represent colour flow in the 
MC "event record". 



Structure Constants of SU(3) 

/l23 = 1 



147 



/246 - 
/l56 



/257 — /: 



345 



/367 — — 



/458 — /f 



678 



2 

V3 
2 



(12) 
(13) 
(14) 

(15) 



Antisymmetric in all indices 
All other fijk = 

(valid for the convention t — ^ ) 
(for the alternative convention t — , multiply all fij by \/2) 



A gluon self-interaction vertex is illustrated in figure 3, to be compared with the quark-gluon 
one in figure 2. We remind the reader that gauge boson self-interactions are a hallmark of non-Abelian 
theories and that their presence leads to some of the main differences between QED and QCD. One 
should also keep in mind that the colour factor for the vertex in figure 3, Ca, is roughly twice as large 
as that for a quark, Cr. 
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Fig. 3: Illustration of a ggg vertex in QCD, before summing/averaging over colours: interaction 
between gluons in the states A^, A'*, and A^ is represented by the structure constant /^^^. 



1.3 The Strong Coupling 

To first approximation, QCD is scale invariant. That is, if one "zooms in" on a QCD jet, one will 
find a repeated self-similar pattern of jets within jets within jets, reminiscent of fractals such as the 
famous Mandelbrot set in mathematics, or the formation of frost crystals in physics. In the context 
of QCD, this property was originally called light-cone scaling, or Bjorken scaling after the famous 
physicist James D. Bjorken. It has since been rebranded by a new generation as conformal invariance, 
a mathematical property of several QCD-"like" theories which are now being studied. It is also closely 
related to the physics of so-called "unparticles", though that is a relation that goes beyond the scope 
of these lectures. 

Regardless of the labeling, if the strong coupling did not run (we shall return to the running of 
the coupling below), Bjorken scaling would be absolutely true. QCD would be a theory with a fixed 
coupling, the same at all scales. This simplified picture already captures some of the most important 
properties of QCD, as we shall discuss presently. 

In the limit of exact Bjorken scaling — QCD at fixed coupling — properties of high-energy in- 
teractions are determined only by dimensionless kinematic quantities, such as scattering angles (pseu- 
dorapidities) and ratios of energy scales'*. For applications of QCD to high-energy collider physics, 
an important consequence of Bjorken scaling is thus that the rate of bremsstrahlung jets with a given 
transverse momentum scales in direct proportion to the hardness of the fundamental partonic scat- 
tering process they are produced in association with. For instance, in the limit of exact scaling, a 
measurement of the rate of 5-GeV jets produced in association with an ordinary Z boson could be 
used as a direct prediction of the rate of 50-GeV jets that would be produced in association with a 
900-GeV Z' boson, and so forth. Our intuition about how many bremsstrahlung jets a given type 
of process is likely to have should therefore be governed first and foremost by the ratios of scales 
that appear in that particular process, as has been highlighted in a number of studies focusing on the 
mass and scales appearing, e.g., in Beyond-the-Standard-Model (BSM) physics processes [6-10]. 
Bjorken scaling is also fundamental to the understanding of jet substructure in QCD, see, e.g., [1 1]. 

In real QCD, the coupling runs logarithmically with the energy, 

2<9as das r., . 



Originally, the observed approximate agreement with this was used as a powerful argument for pointlike substructure in 
hadrons; since measurements at different energies are sensitive to different resolution scales, independence of the absolute 
energy scale is indicative of the absence of other fundamental scales in the problem and hence of pointlike constituents. 
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where the function driving the energy dependence, the beta function, is defined as 



/3{as) = -af (60 + has + 6205 + •••): (17) 
with LO (1-loop) and NLO (2-loop) coefficients 

bo - :^ , (18) 

, UCl-lOTnCAnf-GTRCFUf 153 - 19n/ 

= 24^r^ = 24vr2 " ^^^^ 

Numerically, the value of the strong coupling is usually specified by giving its value at the specific 
reference scale = M^, from which we can obtain its value at any other scale by solving equa- 



tion (16), 



as(Q^) = as(Ml) ^^-^ , (20) 

l + feoa.(M|)lngJ + 0(a2)' 



with relations including the 0{ag) terms available, e.g., in [ ]. Relations between scales not involving 
M| can obviously be obtained by just replacing M| by some other scale Q'"^ everywhere in equa- 
tion (20). As an application, let us prove that the logarithmic running of the coupling implies that an 
intrinsically multi-scale problem can be converted to a single-scale one, up to corrections suppressed 
by two powers of Os, by taking the geometric mean of the scales involved. This follows from expand- 
ing an arbitrary product of individual Os factors around an arbitrary scale /i, using equation (20), 



as(/ii)as(/i2) ■■■a. 



= fl + bo«,lnf ^ \+0{al)) , (21) 

whereby the specific single-scale choice = • • • /^n (the geometric mean) can be seen to push 
the difference between the two sides of the equation one order higher than would be the case for any 
other combination of scales'*. 

The appearance of the number of flavours, Uf, in 60 implies that the slope of the running 
depends on the number of contributing flavours. Since full QCD is best approximated by nj = 3 
below the charm threshold, by nj = 4 from there to the b threshold, and by nj = 5 above that, 
it is therefore important to be aware that the running changes slope across quark flavour thresholds. 
Likewise, it would change across the threshold for top or for any coloured new-physics particles that 
might exist, with a magnitude depending on the particles' colour and spin quantum numbers. 

The negative overall sign of equation (17), combined with the fact that 69 > 0, leads to the 
famous result*' that the QCD coupling effectively decreases with energy, called asymptotic freedom, 
for the discovery of which the Nobel prize in physics was awarded to D. Gross, H. Politzer, and 
F. Wilczek in 2004. An extract of the prize announcement runs as follows: 



""In a fixed-order calculation, the individual scales fit, would correspond, e.g., to the n hardest scales appearing in an 
infrared safe sequential clustering algorithm applied to the given momentum configuration. 

* Perhaps the highest pinnacle of fame for equation (17) was reached when the sign of it "starred" in an episode of the 
TV series "Big Bang Theory". 
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What this year's Laureates discovered was sometliing that, at first sight, seemed com- 
pletely contradictory. The interpretation of their mathematical result was that the 
closer the quarks are to each other, the weaker is the "colour charge". When the 
quarks are really close to each other, the force is so weak that they behave almost 
as free particles". This phenomenon is called "asymptotic freedom". The converse 
is true when the quarks move apart: the force becomes stronger when the distance 
increases'\ 



"More correctly, it is the coupling rather than the force which becomes weak as the distance de- 
creases. The 1 /r^ Coulomb singularity of the force is only dampened, not removed, by the diminishing 
coupling. 

''More correctly, it is the potential which grows, linearly, while the force becomes constant. 

Among the consequences of asymptotic freedom is that perturbation theory becomes better be- 
haved at higher absolute energies, due to the effectively decreasing coupling. Perturbative calculations 
for our 900-GeV Z' boson from before should therefore be slightly faster converging than equivalent 
calculations for the 90-GeV one. Furthermore, since the running of Ug explicitly breaks Bjorken 
scaling, we also expect to see small changes in jet shapes and in jet production ratios as we vary 
the energy. For instance, since high-p^ jets start out with a smaller effective coupling, their intrinsic 
shape (irrespective of boost effects) is somewhat narrower than for \ow-p± jets, an issue which can 
be important for jet calibration. Our current understanding of the running of the QCD coupling is 
summarized by the plot in figure 4, taken from a recent comprehensive review by S. Bethke [ ]. 

As a final remark on asymptotic freedom, note that the decreasing value of the strong coupling 
with energy must eventually cause it to become comparable to the electromagnetic and weak ones, at 
some energy scale. Beyond that point, which may lie at energies of order 10^^ - lO^^GeV (though 
it may be lower if as yet undiscovered particles generate large corrections to the running), we do not 
know what the further evolution of the combined theory will actually look like, or whether it will 
continue to exhibit asymptotic freedom. 

Now consider what happens when we run the coupling in the other direction, towards smaller 
energies. Taken at face value, the numerical value of the coupling diverges rapidly at scales below 1 
GeV, as illustrated by the curves disappearing off the left-hand edge of the plot in figure 4. To make 
this divergence explicit, one can rewrite equation (20) in the following form, 

= — > (22) 
bo In^ 

where 

A ~ 200 GeV (23) 

specifies the energy scale at which the perturbative coupling would nominally become infinite, called 
the Landau pole. (Note, however, that this only parametrizes the purely perturbative result, which 
is not reliable at strong coupling, so equation (22) should not be taken to imply that the physical 
behaviour of full QCD should exhibit a divergence for Q — A.) 

Finally, one should be aware that there is a multitude of different ways of defining both A 
and as{Mz). At the very least, the numerical value one obtains depends both on the renormalization 
scheme used (with the dimensional-regularization-based "modified minimal subtraction" scheme, MS, 
being the most common one) and on the perturbative order of the calculations used to extract them. 
As a rule of thumb, fits to experimental data typically yield smaller values for as{Mz) the higher the 
order of the calculation used to extract it (see, e.g., [12, 13]), with as{Mz)\-LO ^ as(-^^z)|NLO ^ 
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Fig. 4: Illustration of the running of in a theoretical calculation (yellow shaded band) and in 
physical processes at different characteristic scales, from [12]. 



as{Mz)\N'NLO- Further, since the number of flavours changes the slope of the running, the location of 
the Landau pole for fixed as{Mz) depends explicitly on the number of flavours used in the running. 
Thus each value of is associated with its own value of A, with the following matching relations 
across thresholds guaranteeing continuity of the coupling at one loop, 

2_ 2_ 

nf = 4^5 : A5 = A4 (^^^^ A4 = A5 (1^)'' , (24) 

2_ 2_ 

(A q \ 25 / TO \ 27 

^'^^^[lu) ■ ^^^^ 

It is sometimes stated that QCD only has a single free parameter, the strong coupling. Appeal- 
ing as this may be, it is a bit of an overstatement. Even in the perturbative region, the beta function 
depends explicitly on the number of quark flavours, as we have seen, and thereby also on the quark 
masses. Furthermore, in the non-perturbative region around or below Aqcd> the value of the perturba- 
tive coupling, as obtained, e.g., from equation (22), gives little or no insight into the behaviour of the 
full theory. Instead, universal functions (such as parton densities, form factors, fragmentation func- 
tions, etc), effective theories (such as the Operator Product Expansion, Chiral Perturbation Theory, or 
Heavy Quark Effective Theory), or phenomenological models (such as Regge Theory or the String 
and Cluster Hadronization Models) must be used, which in turn depend on additional non-perturbative 
parameters whose relation to, e.g., as{Mz), is not a priori known. For some of these questions, such 
as hadron masses, lattice QCD can furnish important additional insight, but for multi-scale and/or 
time-evolution problems, the applicability of lattice methods is still severely restricted. 
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Fig. 5: Left: Rutherford scattering of quarks in QCD, exemplifying the type of process that dominates 
the short-distance interaction cross section at hadron colliders. Right: an example of what such a 
reaction may look like in a detector, in this case the ATLAS experiment. 



2 Perturbative QCD 

Our main tool for solving QCD for high-energy collider physics is perturbative quantum field theory, 
the starting point for which is Matrix Elements (MEs) which can be calculated systematically at fixed 
orders in the strong coupling Og- At least at lowest order (LO), the procedure is standard textbook 
material [5] and it has also by now been highly automated, by the advent of tools like CalcHep [14], 
CompHep [ ], MadGraph [16], and others [17-21]. Here, we require only that the reader has a 
basic familiarity with the methods involved from graduate-level particle physics courses based, e.g., 
on [5,21]. Our main concern are the uses to which these calculations are put, their limitations, and 
ways to improve on the results obtained with them. 

For illustration, take one of the most commonly occurring processes in hadron collisions — 
Rutherford scattering of two quarks via a t-channel gluon exchange — which has the differential 
cross section 

, , da TT 4 2 s'^ + u'^ 
m^m : ^ = p9«.^^, (26) 

with the 2 — 2 Mandelstam variables ("hatted" to emphasize that they refer to a partonic 2 — 2 
scattering rather than the full pp — jets process) 

s = {Pi+P2f, (27) 

,9 . (1 — COS 6) 
t = {P3-pif = -s^- ^ -, (28) 

u = {pA-piY = -s^ 2 '- . (29) 

This process is illustrated in the left-hand pane of figure 5, including a rough (formally leading- A'^c') 
representation of the "colour transfer" mediated by the gluon (as was discussed in Section 1.2). 

Reality, however, is more complicated; the picture on the right-hand pane of figure 5 shows 
a real dijet event, as recorded by the ATLAS experiment. The complications to be addressed when 
going from left to right in figure 5 are: firstly, additional jets, a.k.a. real-emission corrections, which 
significantly change the topology of the final state, potentially shifting jets in or out of an experimen- 
tally defined acceptance region. Secondly, loop factors, a.k.a. virtual corrections, change the number 
of available quantum paths through phase space, and hence modify the normalization of the cross 
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section (total and differential). And finally, additional corrections to the simple factorized pertur- 
bative picture are generated by components such as hadronization and the underlying event. These 
corrections must be taken into account to complete our understanding of QCD and connect the short- 
distance physics with macroscopic experiments. Apart from the perturbative expansion itself, the 
most powerful tool we have to organize this vast calculation, is factorization. 

2.1 Factorization 

When applicable, factorization allows us to subdivide the calculation of an observable into a pertur- 
batively calculable short-distance part and an approximately universal long-distance part, the latter 
of which may be modeled and constrained by fits to data. Factorization can also be applied multiple 
times, to break up a complicated calculation into simpler pieces that can be treated as approximately 
independent, such as when dealing with successive emissions in a parton shower, or when factoring 
off decays of long-lived particles from a hard production process. 

Using coUinear factorization (see, e.g., [4,23]), the differential cross section for an observable 
O in hadron-hadron collisions can be computed as: 

^ = Y1 /'dx.dxfe ^ / dcl>^ ft{xa,^^F)f'S\x,,^^F) '^^^DpiO ^ (30) 

a,b h p J dO 

where the outer sum runs over all partonic constituents, a and h of the colliding hadrons, /ii 2, re- 
spectively, and the inner sum runs over all possible final states, ah ^ F (with the standard final-state 
phase-space differential denoted d^^ ). 

Before we discuss the integrand — composed of the factors da , and Dp — let us first re- 
emphasize the crucial feature of equation (30); it separates the calculation of the cross section into two 
independent pieces, one of which is the perturbatively calculable short-distance cross section, da , and 
the other of which is the product of parton distribution functions (PDFs), fafb, with a fragmentation 
function (FF), Dp, with the latter components being universal functions'' whose forms are a priori 
unknown but which can be constrained in one process and then reused in another. The dividing line 
between the two is drawn at an arbitrary ("user-defined") scale jip, called the. factorization scale. 

Returning now to the integrand, the parton density functions, fj'{xj,^F), parametrize the 
distribution of partons of type j carrying momentum fraction Xj inside a hadron of type hi when 
probing the latter at the factorization scale fip. (Note: issues specific to PDFs in the context of Monte 
Carlo event generators will be covered in Section 3.1.) The partonic scattering cross section daab-^F 
is calculable in fixed-order perturbation theory as 

daab^F =:r^\Mab~^F\'^{^p;iJ'F,fJ'R), , (31) 

with the matrix element squared for the process ab — )• F, appropriately summed and averaged 
over helicities and/or colours, and evaluated at the factorization and renormalization scales /x p and fin, 
respectively. The fragmentation functions (FFs), Dp{0 — )• 0,fip) parametrize the transition from 
partonic final state to the hadronic observable (bremsstrahlung, hadronization, jet definition, etc). 

'At least, they are universal within the framework of collinear factorization. In full QCD, there are several types of 
corrections, including also some perturbative ones, that go beyond this framework, such as small-a; effects and multiple 
parton interactions, both of which mandate the introduction of objects that go beyond the scope of collinear-factorized 
PDFs. In the case of small-a:: evolution, these more general objects are so-called unintegrated PDFs, which have an explicit 
dependence on the parton transverse momentum in addition to the factorization scale, while multi-parton interactions require 
explicit multi-parton and/or generalized (impact-parameter-dependent) PDFs. 
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There is some arbitrariness involved in this division of the calculation into a short-distance 
and a long-distance part. Firstly, one has to choose a value for the dividing scale, ^p. Some heuristic 
arguments to guide in the choice of factorization scale are the following. On the long-distance side, the 
PDFs include a (re)summation of multiple emissions (bremsstrahlung) all the way up to the scale fip- 
It would therefore not make much sense to take fi p significantly larger than the scales characterizing 
resolved particles on the short-distance side of the calculation (i.e., the particles appearing explicitly 
in ^p); otherwise the PDFs would be including sums over radiations as hard as or harder than those 
included explicitly in the matrix element which would result in double-counting. On the other hand, 
it should not be taken much lower than the scales appearing in the matrix element either, since, as 
we shall see in subsequent chapters, fixed-order matrix elements are at most able to include part of 
such multiple-bremsstrahlung emissions, and hence a low choice of factorization scale would lead to 
problems with "undercounting" of such corrections. 

For matrix elements characterized by a single well-defined scale, such as the scale in deeply 
inelastic scattering (DIS) or the invariant-mass scale s in Drell-Yan production {qq — )• Z/7* — )• ), 
such arguments essentially fix the preferred scale choice, which may then be varied by a factor of 2 
(or larger) around the nominal value in order to estimate uncertainties. For multi-scale problems, 
however, such as pp — )• Z/W + n jets, there are several a priori equally good choices available, from 
the lowest to the highest QCD scales that can be constructed from the final-state momenta, usually 
with several dissenting groups of theorists arguing over which particular choice is best. Suggesting 
that one might simply measure the scale would not really be an improvement, as the factorization scale 
is fundamentally unphysical and therefore unobservable (similarly to gauge or convention choices). 
One plausible strategy is to look at higher-order (NLO or NNLO) calculations, in which correction 
terms appear that explicitly remove the over- or under-counting introduced by the initial scale choice 
up to the given order, thus reducing the overall dependence on it and stabilizing the final result. From 
such comparisons, a "most stable" initial scale choice can in principle be determined, which then 
furnishes a reasonable starting point, but we emphasize that the question is intrinsically ambiguous, 
and no "golden recipe" is likely to magically give all the right answers. The best we can do is to vary 
the value of fip not only by an overall factor, but also by exploring different possible forms for its 
functional dependence on the momenta appearing in ^p. In this way, one could hope to provide a 
more complete uncertainty estimate for multi-scale problems. 

Secondly, and more technically, at NLO and beyond one also has to settle on a factorization 
scheme in which to do the calculations. For all practical purposes, students focusing on LHC physics 
are only likely to encounter one such scheme, the modified minimal subtraction (MS) one already 
mentioned in the discussion of the definition of the strong coupling in Section 1.3. At the level of 
these lectures, we shall therefore not elaborate further on this choice here. 

2.2 Infrared Safety 

The second perturbative tool, infrared safety, provides us with a special class of observables which 
have minimal sensitivity to long-distance physics, and which can be consistently computed in per- 
turbative QCD (pQCD). By "infrared", we here mean any limit that involves a low scale (i.e., any 
non-UV limit), without regard to whether it is collinear or soft*^. An observable is infrared safe if: 

1 . (Safety against soft radiation): Adding any number of infinitely soft particles should not change 
the value of the observable. 

^This distinction will be discussed further in Section 2.4. 



14 



2. (Safety against collinear radiation): Splitting an existing particle up into two comoving parti- 
cles, with arbitrary fractions z and I — z, respectively, of the original momentum, should not 
change the value of the observable. 



If both of these conditions are satisfied, any long-distance non-perturbative corrections will be sup- 
pressed by the ratio of the long-distance scale to the short-distance one to some (observable-dependent) 
power, typically 



where Quv denotes a generic hard scale in the problem, and Qir ~ Aqcd ~ 0(1 GeV). 

Due to this power suppression, IR safe observables are not so sensitive to our lack of ability to 
solve the strongly coupled IR physics, unless of course we go to processes for which the relevant hard 
scale, QuY, is small (such as minimum-bias, soft jets, or small-scale jet substructure). Even when 
a high scale is present, however, as in resonance decays, jet fragmentation, or underlying-event-type 
studies, infrared safety only guarantees us that infrared corrections are small, not that they are zero. 
Thus, ultimately, we run into a precision barrier even for IR safe observables, which only a reliable 
understanding of the long-distance physics itself can address. 

To constrain models of long-distance physics, one needs infrared sensitive observables''. In- 
stead of the suppressed corrections above, the perturbative prediction for such observables contains 
logarithms 



which grow increasingly large as Qm/Q\jv — ^ 0. As an example, consider such a fundamental 
quantity as particle multiplicities; in the absence of nontrivial infrared effects, the number of partons 
that would be mapped to hadrons in a naive local-parton-hadron-duality [24] picture would tend log- 
arithmically to infinity as the IR cutoff is lowered. Similarly, the distinction between a charged and 
a neutral pion only occurs in the very last phase of hadronization, and hence observables that only 
include charged tracks, for instance, are always IR sensitive'*'. 

Two important categories of infrared safe observables that are widely used are event shapes and 
jet algorithms. Jet algorithms are perhaps nowhere as pedagogically described as in last year's ESHEP 
lectures by Salam [ , Chapter 5]. Event shapes in the context of hadron colliders have not yet been as 
widely explored, but the basic phenomenology is introduced also by Salam and collaborators in [ ], 
with a first measurement reported by CMS [27] and a proposal to use them also for the characterization 
of minimum-bias events put forth in [28]. 

Let us here merely emphasize that the real reason to prefer infrared safe jet algorithms over 
unsafe ones is not that they necessarily give very different or "better" answers in the experiment 
— experiments are infrared safe by definition, and the difference between infrared safe and unsafe 
algorithms may not even be visible when running the algorithm on experimental data — but that it is 
only possible to compute perturbative QCD predictions for the infrared safe ones. Any measurement 
performed with an infrared unsafe algorithm can only be compared to calculations that include a 
detailed hadronization model. This both limits the number of calculations that can be compared to 

' Hence it is not always the case that infrared safe observables are preferable — the purpose decides the tool. 
'"This remains true in principle even if the tracks are clustered into jets, although the energy clustered in this way does 
provide a lower bound on Quv in the given event, since "charged + neutral > charged-only". 



IR Safe Observables: IR corrections oc 




(32) 



IR Sensitive Observables: IR Corrections oc q" log' 




m <2n 



(33) 



15 



and also adds an a priori unknown sensitivity to the details of the hadronization description, details 
which one would rather investigate and constrain separately, in the framework of more dedicated 
fragmentation studies. 

2.3 Fixed-Order QCD: Matrix Elements 

Schematically, we express the all-orders differential cross section for an observable O, in the pro- 
duction of F -I- anything (= inclusive F production, with F an arbitrary final state), in the following 
way: 



where, for compactness, we have suppressed all PDF and luminosity normalization factors. The sum 
over k represents a sum over additional "real-emission" corrections, called legs, and the sum over i 
runs over additional virtual corrections, loops. Without the 5 function, the formula would give the total 
integrated cross section, instead of the cross section differentially in O. The purpose of the 5 function 
is thus to project out hypersurfaces of constant O in the full d<I>^^^ phase space, with C'(<^>^^^) a 
function that defines O evaluated on each specific momentum configuration, ^p_^^f.- 

We recover the various fixed-order truncations of pQCD by limiting the nested sums in equa- 
tion (34) to include only specific values of k + i. Thus, 

k = 0, £ = =^ Leading Order (usually tree-level) for inclusive F production 
k = n, £ = =^ Leading Order for F + n jets 

k + £<n, =^ N"LO for F (includes N"-iLO for F + Ijet, N"-2lo for F + 



Already at this stage, before entering into the details of the calculations, we can make several obser- 
vations on how numerical values of cross sections and decay widths must be computed in fixed-order 
perturbation theory. 

Firstly, the dimensionality of the phase space to be integrated increases by d = 3 for each leg 
we add. In dimensions higher than 5, the fastest converging numerical integration algorithm is Monte 
Carlo integration [29], whose purely stochastic error oc 0{1/VJ7), with J\f the number of generated 
points, is independent of dimension, while all other algorithms scale with powers of the dimension. 
Therefore, virtually all numerical cross section calculations are based on Monte Carlo techniques in 
one form or another, the simplest being the Rambo algorithm [30] which can be expressed in about 
half a page of code and generates a fiat scan over n-body phase space' ' . 

Secondly, due to the infrared singularities in perturbative QCD, the functions to be integrated, 
|A^p, are highly non-uniform for large k, which implies that we will have to be clever in the way 
we sample phase space if we want the integration to converge in any reasonable amount of time — 
simple algorithms like Rambo quickly become inefficient for k greater than a few. To address this 
bottleneck, the simplest step up from RAMBO is to introduce generic (i.e., automated) importance- 
sampling methods, such as offered by the VEGAS algorithm [ , ]. This is still the dominant basic 
technique, although most modern codes do employ several additional refinements, such as several 
different copies of VEGAS running in parallel (multi-channel integration), to further optimize the 

"strictly speaking, RAMBO is only truly uniform for massless particles. Its massive variant makes up for phase-space 
biases by returning weighted momentum configurations. 




(34) 



s legs 



s loops 



2jets, and so on up to LO for F + njets) . 
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sampling. Alternatively, a few algorithms incorporate the singularity structure of QCD explicitly 
in their phase-space sampling, either by directly generating momenta distributed according to the 
leading-order QCD singularities, in a sort of "QCD-preweighted" analog of Rambo, called Sarge 
[33], or by using all-orders Markovian parton showers to generate them (ViNCiA [ , . ]). 

Thirdly, for A; > 1,£ = 0, we are really not considering inclusive F production anymore; 
instead, we are considering the LO contribution to the process F + A; jets. However, if we simply 
integrate over all momenta, as implied by the integration over d^p_^_f^ in equation (34), we would 
be including configurations in which one or more of the k partons become coUinear or soft, leading 
to singularities in the integration region. At the LO level, this problem can only be mitigated by re- 
stricting the integration region to only include "hard, well-separated" momenta. As discussed above, 
due to the approximate Bjorken scaling of QCD, it would be meaningless to express this requirement 
in dimensionful terms, as an absolute scale. Instead, it is the ratios of scales present in any given 
process that determine whether such enhancements are present or absent: a 50-GeV jet would be con- 
sidered hard and well-separated if produced in association with an ordinary Z boson, while it would 
be considered soft if produced in association with a 900-GeV Z' boson [6-8]. Thus, for example, it 
would be a complete disaster to use the same dimensionful phase-space cuts for Z'+jets as one uses 
for Z+jets (unless of course the Z' happens to have a mass scale very close to the Z one). A good 
rule of thumb is that if cjfc+i ~ Ufc (at whatever order you are calculating), then you are integrating 
over a region in which the perturbative series is no longer converging, or is converging too slowly for 
a fixed-order truncation of it to be reliable. For fixed-order perturbation theory to be applicable, you 
must have (Tk+i ^ afc. In the discussion of parton showers and resummations in Section 2.4, we shall 
see how the region of applicability of perturbation theory can be extended. 

And finally, the virtual amplitudes, for £ > 1, are divergent for any point in phase space. 
However, as encapsulated by the famous KLN theorem [ , ], unitarity (which essentially expresses 
probability conservation) puts a powerful constraint on the IR divergences'", forcing them to cancel 
exactly against those coming from the unresolved emissions that we had to cut out above, order by 
order, making the complete answer for fixed k + £ = n finite. Nonetheless, since this cancellation 
happens between contributions that formally live in different phase spaces, a main aspect of loop- 
level higher-order calculations is how to arrange for this cancellation in practice, either analytically or 
numerically, with many different methods currently on the market. 

A convenient way of illustrating the terms of the perturbative series that a given matrix-element- 
based calculation includes is given in figure 6. In the left-hand pane, the shaded box corresponds to 
the lowest-order "Born-level"'^ matrix element squared. This coefficient is non-singular and hence 
can be integrated over all of phase space, which we illustrate by letting the shaded area fill all of 
the relevant box. A different kind of leading-order calculation is illustrated in the right-hand pane of 
figure 6, where the shaded box corresponds to the lowest-order matrix element squared for F + 2 jets. 
This coefficient diverges in the part of phase space where one or both of the jets are unresolved (i.e., 
soft or coUinear), and hence integrations can only cover the hard part of phase space, which we reflect 
by only shading the upper half of the relevant box. 

Figure 7 illustrates the inclusion of NLO virtual corrections. To prevent confusion, first a point 
on notation: by (7q'\ we intend 




(35) 
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The loop integrals also exhibit UV divergences, but these are deah with by renormalization. 
Photo from nobelprize . org 
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Fig. 6: Coefficients of tlie perturbative series covered by LO calculations. Left: F production at 
lowest order. Right: F + 2 jets at LO, with the half-shaded box illustrating the restriction to the region 
of phase space with exactly 2 resolved jets. The total power of as for each coefficient is n = k + £. 
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Fig. 7: Coefficients of the perturbative series covered by NLO calculations. Left: F production at 
NLO. Right: F + Ijet at NLO, with half-shaded boxes illustrating the restriction to the region of 
phase space with exactly 1 resolved jet. The total power of for each coefficient is n = k + £. 



which is of order as relative to the Born level. Compare, e.g., with the expansion of equation (34) to 
order k + i = 1. In particular, should not be confused with the integral over the 1-loop matrix 
element squared (which would be of relative order and hence forms part of the NNLO coefficient 

(2) 

(Jq ). Returning to figure 7, the unitary cancellations between real and virtual singularities imply that 
we can now extend the integration of the real correction in the left-hand pane over all of phase space, 
while retaining a finite total cross section, 

J J J (36) 

where the divergence caused by integrating the third term over all of phase space is canceled by that 
coming from the integration over loop momenta in the second term. However, if our starting point for 
the NLO calculation is a process which already has a non-zero number of hard jets, we must continue 
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Fig. 8: Coefficients of tlie perturbative series covered by an NNLO calculation. The total power of 
as for each coefficient is n = k + i. Green shading represents the full perturbative coefficient at the 
respective k and i. 



to impose that at least that number of jets must still be resolved in the final-state integrations, 

ar°(p±min) = / d<l>, \Mf^\' + [ d$i 2Re[A^«A^f)*] + / d^^ l-M^ I' 

^PX>P±min •^P±>PXmin •^PXl>P±min 

= Crf'\p± > p±nnn) + (^[^\p± > P±mm) + (^i^\p±l > P±mm) , 

(37) 

where the restriction to at least one jet having p± > p±mm has been illustrated in the right-hand pane 
of figure 7 by shading only the upper part of the relevant boxes. In the last term in equation (37), the 
notation p±i is used to denote that the integral runs over the phase space in which at least one "jet" 
(which may consist of one or two partons) must be resolved with respect to p±ram- Here, therefore, an 
explicit dependence on the algorithm used to define "a jet" enters for the first time. This is discussed 
in more details in the ESHEP lectures by Salam [ ]. 

To extend the integration to cover also the case of 2 unresolved jets, we must combine the left- 
and right-hand parts of figure 7 and add the new coefficient 

= |A^«P + 2Re[>l(^Wr], (38) 
as illustrated by the diagram in figure 8. 

2.4 Infinite-Order QCD: Parton Showers 

In the preceding section, we noted two conditions that had to be valid for fixed-order truncations of 
the perturbative series to be valid: firstly, the strong coupling must be small for perturbation theory 
to be valid at all. This restricts us to the region in which all scales Qi » Aqcd- We shall maintain 
this restriction in this section, i.e., we are still considering perturbative QCD. Secondly, however, in 
order to be allowed to truncate the perturbative series, we had to require cjfc+i ^ au, i.e., the correc- 
tions at successive orders must become successively smaller, which — due to the enhancements from 
soft/collinear singular (conformal) dynamics — effectively restricted us to consider only the phase- 
space region in which all jets are "hard and well-separated", equivalent to requiring all Qi/Qj w 1. 
In this section, we shall see how to lift this restriction, extending the applicability of perturbation 
theory into regions that include scale hierarchies, Qi » Qj ^> Aqcd> such as occur for soft jets, jet 
substructure, etc. 
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In fact, the simultaneous restriction to all resolved scales being larger than Aqcd cind no large 
hierarchies is extremely severe, if taken at face value. Since we collide and observe hadrons (— )• low 
scales) while simultaneously wishing to study short-distance physics processes (— )• high scales), it 
would appear trivial to conclude that fixed-order pQCD is not applicable to collider physics at all. So 
why do we still use it? 

The answer lies in the fact that we actually never truly perform a fixed-order calculation in 
QCD. Let us repeat the factorized formula for the cross section, equation (30), 

% = Y. I'dXadXk [ dcD^ f^'{Xa,^^F)f^%Xk,^^F) ^^^1)^(0 ^ O, ^f) • (39) 
a,b Jo p J dO 

Although daab^F does represent a fixed-order calculation, the parton densities, /^^ and /^^, include 
so-called resummations of perturbative corrections to all orders from the initial scale of order the 
mass of the proton, up to the factorization scale, fip- Note that the oft-stated mantra that the PDFs are 
purely non-perturbative functions is therefore misleading. True, they are defined as essentially non- 
perturbative functions at some very low scale, but, if fip is taken large, they necessarily incorporate 
a significant amount of perturbative physics as well. On the "fixed-order side", all we have left to 
ensure in daab^F is then that there are no large hierarchies remaining between fip and the QCD 
scales appearing in ^p. Likewise, in the final state, the fragmentation functions. Dp, include infinite- 
order resummations of perturbative corrections all the way from fi p down to some low scale, with 
similar caveats concerning mantras about their non-perturbative nature as for the PDFs. 



2.4.1 Step One: Infinite Legs 

The infinite-order resummations that are included in objects such as the PDFs and FFs in equation (39) 
(and in their parton-shower equivalents) rely on some very simple and powerful properties of gauge 
field theories. One way to arrive at them is the following; assume we have computed the Born-level 
cross section for some process, F, and that this process contains some number of coloured partons'"^. 
For each pair of (massless) colour-anticolour charges A and B in F, it is then a universal property of 
QCD that the cross sections for F + 1 partons, da^p]^^ will include a factor 

dafl^ = gl i MAB-^aib'^^^'^^^ + less singular terms ) x da^ , (40) 
^ V Sal Sl6 / 

where, for compactness, we have lumped some uninteresting normalization factors''' into MAB^aib, 
gl = 47ras is the strong coupling, a and h represent partons A and B after the branching (i.e., they 
include possible recoil effects) and sn is the invariant between parton i and the emitted "-i-l" parton. 
Intuitively, this structure follows from the simple observations illustrated by the left and middle panes 
of figure 9; the Feynman diagram in which parton "1" is emitted from the "o" (or "6") leg has a pole 
for Sal — ^ (sift — )• 0), corresponding to the intermediate propagator "a*" ("6*") going on shell 
(middle pane). Summing the two and squaring them, i.e., including their mutual interference, one 
obtains the structure in equation (40), which is called the Eikonal factor. 

Assume further that octet colour charges (gluons) may be represented as the sum of a colour triplet and an antitriplet 
charge — compare, e.g., with the illustrations of gluon colour flow, Figures 2 and 3. This picture of octets is correct up to 
corrections of order 1/Nq, which will be good enough for our purposes here. 

" I.e., JVAB^aib contains colour and phase-space normalization factors. Up to mildly non-universal corrections of order 
1/Nq (which depend on whether the emitting particles are quarks or gluons), it is A/'As->ai6 = 2CA/(167r^) 
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a) Original Configuration: dap^ b) A contribution to dap!^_^ c) Recursion — ^ dcr^J^j,^ 

Fig. 9: a) and b) Illustration of the QCD singularities induced by on-shell propagators, c) The approx- 
imation obtained in the first step can be iterated to add additional legs. 



.(0) 



The leading part of the singularity structure to which we have already referred many times is 
clearly visible here: if we integrate over the entire phase space including the region Sai — s- 0, — > 0, 
we end up with a double pole. If we instead regulate the divergence by cutting off the integration at 
some minimal perturbative cutoff scale nfj^, we end up with a logarithm squared of that scale'^. This 
is a typical example of "large logs" being generated by the presence of scale hierarchies. 

Before we continue, it is worth noting that equation (40) is often rewritten in other forms to 
emphasize specific aspects of it. One such rewriting is thus to reformulate the invariants sn appearing 
in equation (40) in terms of energies and angles, 

Sij = 2EiEj (1 - cos%) . (41) 

Rewritten in this way, the differentials in equation (40) become 

dsai dsit, dEi d0ai , d£'i d^if, 
°^ ~^ a ^ T' ^ ■ 

■Sal Sib ^1 C'al -C/1 fc'lb 

This kind of rewriting enables an intuitively appealing categorization of the singularities as related to 
vanishing energies and angles, called soft and coUinear limits, respectively. Although such formu- 
lations have undeniably been helpful in obtaining many important results in QCD, one should still 
keep in mind that Lorentz non-invariant formulations come with similar caveats and warnings as do 
gauge non-invariant formulations of quantum field theory: while they can be practical to work with 
at intermediate stages of a calculation, one should be careful with any physical conclusions that rely 
explicitly on them. We shall therefore here restrict ourselves to a Lorentz invariant formalism based 
directly on equation (40). The collinear limit is then replaced by a more general single-pole limit 
in which a single parton-parton invariant vanishes (as, for instance, when a pair of partons become 
collinear), while the soft limit is replaced by one in which two (or more) invariants involving the same 
parton vanish simultaneously (as, for instance by that parton becoming soft in a frame defined by two 
or more hard partons). This avoids frame-dependent ambiguities from entering into the language, at 
the price of a slight reinterpretation of what is meant by collinear and soft. 

Independently of rewritings and philosophy, the real power of equation (40) lies in the fact that 
it is universal. Thus, for any process F, we can apply equation (40) in order to get an approximation 
for dcjp^^i . We may then, for instance, take our newly obtained expression for F + 1 as our arbitrary 
process and crank equation (40) again, to obtain an approximation for dap +2 , and so forth. What we 
have here is therefore a very simple recursion relation that can be used to generate approximations 

"The precise definition of fif-i^ is not unique. Any scale choice that properly isolates the singularities from the rest of 
phase space will do, with some typical choices being, for example, invariant-mass and/or transverse-momentum scales. 
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to leading-order cross sections with arbitrary numbers of additional legs. The quality of this approx- 
imation is governed by how many terms besides the leading one shown in equation (40) are included 
in the game. Including all possible terms, the most general form for the cross section at F + n jets, 
restricted to the phase-space region above some infrared cutoff scale ;Uir, has the following algebraic 
structure, 

aj^l^ = (ln2" + In^""! + ln2"-2 + . . . + In + 7^) (43) 

where we use the notation In'^ without an argument to denote generic functions of transcendentality A 
(the logarithmic function to the power A being a "typical" example of a function with transcendentality 
A appearing in cross section expressions, but also dilogarithms and higher logarithmic functions'^ of 
transcendentality > 1 should be implicitly understood to belong to our notation In^). The last term, IZ, 
represents a rational function of transcendentality 0. We shall also use the nomenclature singular and 
finite for the In'^ and Ti terms, respectively, a terminology which reflects their respective behaviour in 
the limit jiin — > 0. 

The simplest approximation one can build on equation (43), dropping all but the leading In^*^ 
term in the parenthesis, is thus the leading-transcendentality approximation. This approximation is 
better known as the DLA (double logarithmic approximation), since it generates the correct coeffi- 
cient for terms which have two powers of logarithms for each power of as, while terms of lower 
transcendentalities are not guaranteed to have the correct coefficients. In so-called LL (leading- 
logarithmic) parton shower algorithms, one generally expects to reproduce the correct coefficients 
for the In^*^ and In^"^'^ terms. In addition, several formally subleading improvements are normally 
also introduced in such algorithms (such as explicit momentum conservation, gluon polarization and 
other spin-correlation effects [ ], higher-order coherence effects, renormalization scale choices [ ], 
finite- width effects [ ], etc), as a means to improve the agreement with some of the more subleading 
coefficients as well, if not in every phase-space point then at least on average. Though LL showers do 
not magically acquire NLL (next-to-leading-log) precision from such procedures, one therefore still 
expects a significantly better average performance from them than from corresponding "strict" LL 
analytical resummations. A side effect of this is that it is often possible to "tune" shower algorithms 
to give better-than-nominal agreement with experimental distributions, by adjusting the parameters 
controlling the treatment of subleading effects. One should remember, however, that there is a limit 
to how much can be accomplished in this way — at some point, agreement with one process will only 
come at the price of disagreement with another, and at this point further tuning would be meaningless. 

Applying such an iterative process on a Bom-level cross section, one obtains the description of 
the full perturbative series illustrated in the left-hand pane of figure 10. The yellow (lighter) shades 
used here for A; > 1 indicate that the coefficient obtained is not the exact one, but rather an approxima- 
tion to it that only gets its leading singularities right. However, since this is still only an approximation 
to infinite-order tree-level cross sections (we have not yet included any virtual corrections), we cannot 
yet integrate this approximation over all of phase space, as illustrated by the yellow boxes being only 
half filled on the left-hand side of figure 10; the summed total cross section would still be infinite. 
This particular approximation would therefore still appear to be very useless indeed — on one hand, 
it is only guaranteed to get the singular terms right, but on the other, it does not actually allow us to 
integrate over the singular region. In order to obtain a truly all-orders calculation, the constraint of 
unitarity must also be explicitly imposed, which furnishes an approximation to all-orders loop correc- 
tions as well. Let us therefore emphasize that figure 10 is included for pedagogical purposes only; all 

"Note: due to the theorems that allow us, for instance, to rewrite dilogarithms in different ways with logarithmic and 
lower "spillover" terms, the coefficients at each A are only well-defined up to reparametrization ambiguities involving the 
terms with transcendentality greater than A. 
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Fig. 10: Coefficients of tlie perturbative series covered by LO + LL calculations, Left: without impos- 
ing unitarity, and Right: imposing unitarity order by order for each n = k + i. Green (darker) shading 
represents the full perturbative coefficient at the respective k and £. Yellow (lighter) shading repre- 
sents an LL approximation to it. Half-shaded boxes indicate phase spaces in which we are prohibited 
from integrating over the IR singular region, as discussed in sections 2.3 and 4.2. 



resummation calculations, whether analytical or parton-shower based, include virtual corrections as 
well and consequently yield finite total cross sections, as will now be described. 

2.4.2 Step Two: Infinite Loops 

Order-by-order unitarity, such as used in the KLN theorem, implies that the singularities caused by 
integration over unresolved radiation in the tree-level matrix elements must be canceled, order by 
order, by equal but opposite-sign singularities in the virtual corrections at the same order. That is, 
from equation (40), we immediately know that the 1-loop correction to dap must contain a term, 

dcxL^^ = -gsMAB^aibdaP f '^^(^^ _^ less singular terms, (44) 

J Sal Sib 

that cancels the divergence coming from equation (40) itself. Further, since this is universally true, we 
may apply equation (44) again to get an approximation to the corrections generated by equation (40) 
at the next order and so on. By adding such terms explicitly, order by order, we may now bootstrap our 
way around the entire perturbative series, using equation (40) to move horizontally and equation (44) 
to move along diagonals of constant n = k + £. Since real-virtual cancellations are now explicitly 
restored, we may finally extend the integrations over all of phase space, resulting in the picture shown 
on the right-hand pane of figure 10. 

The right-hand pane, not the left-hand one, corresponds to what is actually done in resummation 
calculations, both of the analytic and parton-shower types Physically, there is a significant and 
intuitive meaning to the imposition of unitarity, as follows. 

Take a jet algorithm, with some measure of jet resolution, Q, and apply it to an arbitrary sample 
of events, say dijets. At a very crude resolution scale, corresponding to a high value for Q, you find 
that everything is clustered back to a dijet configuration, and the 2-jet cross section is equal to the total 
inclusive cross section, 

O"tot = 0'F;incl • (45) 

'**In the way these calculations are formulated in practice, they in fact rely on one additional property, called exponenti- 
ation, that allows us to move along straight vertical lines in the loops-and-legs diagrams. However, since the two different 
directions furnished by equations (40) and (44) are already sufficient to move freely in the full 2D coefficient space, we 
shall use exponentiation without extensively justifying it here. 
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At finer resolutions, decreasing Q, you see tliat some events tliat were previously classified as 2-jet 
events contain additional, lower-scale jets, that you can now resolve, and hence those events now 
migrate to the 3-jet bin, while the total inclusive cross section of course remains unchanged. 



ftot = <7F;mcl = 0"F;excl(Q) + 0'F+l;incl(Q) 



(46) 



where "incl" and "excl" stands for inclusive and exclusive cross sections'^, respectively, and the Q- 
dependence in the two terms on the right-hand side must cancel so that the total inclusive cross section 
is independent of Q. Later, some 3-jet events now migrate further, to 4 and higher jets, while still more 
2-jet events migrate into the 3-jet bin, etc. For arbitrary n and Q, we have 



This equation expresses the trivial fact that the cross section for n or more jets can be computed as 
the total inclusive cross section for F minus a sum over the cross sections for F + exactly m jets 
including all m < n. On the theoretical side, it is these negative terms which must be included in the 
calculation, for each order n = k + i,to restore unitarity. Physically, they express that, at a given scale 
Q, a given event will be classified as having either 0, 1, 2, or whatever jets. Or, equivalently, for each 
event we gain in the 3-jet bin as Q is lowered, we must loose one event in the 2-jet one; the negative 
contribution to the 2-jet bin is exactly minus the integral of the positive contribution to the 3-jet one, 
and so on. We may perceive of this detailed balance as an evolution of the event structure with Q, for 
each event, which is effectively what is done in parton-shower algorithms, to which we shall return in 
Section 4.1. 



In a complete high-energy collision, many different physics (sub-)processes contribute to the total 
observed activity. We here give a very brief overview of the main aspects of soft QCD that are relevant 
for hadron-hadron collisions, such as parton distribution functions, minimum-bias and soft-inclusive 
physics, and the so-called "underlying event". This will be kept at a strictly pedestrian level and is 
largely based on the review in [41]. A discussion of the modeling of these components, as well as 
a discussion of the process of hadronization, is deferred to the relevant parts of Section 4 on Monte 
Carlo event generators. 

3.1 Parton Densities 

Physically, parton densities express the fact that hadrons are composite, with a time-dependent struc- 
ture, illustrated in figure 11. More formally, they are defined by the factorization theorem discussed 
in Section 2.1. Occasionally, the words structure functions and parton densities are used interchange- 
ably. However, there is a very important distinction between the two, which we find often in (quantum) 
physics: one is a physical observable, the other is a "fundamental" quantity extracted from it. 

Structure functions, such as F2, are completely unambiguous physical observables, which can 
be measured, for instance, in DIS processes. (For a definition, see, e.g., [ ].) From these, and 
other observables, a set of more fundamental and theoretically useful objects, parton density func- 
tions (PDFs), can be extracted, but there is a price; since the parton densities are not, themselves, 

^^F inclusive = F plus anything. F exclusive = F and only F. Thus, (TF;inci = Jlli^o ""f'+fcicxci 




(47) 
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Fig. 1 1 : Illustration (from [ ]) of partonic fluctuations inside a proton beam. 



physically observable, they can only be defined within a specific factorization scheme, order by or- 
der in perturbation theory. The only exception is at leading order, at which they have a very simple 
physical interpretation, as the probability of finding a quark of a given flavor and carrying a given 
momentum fraction, x, inside a hadron of a given type, probed at a specific scale, Q^. They are then 
related to the structure function F2 by their charge-weighted momentum sum. 



where fi denotes the parton density for a parton of flavor/type i. When going to higher orders, we 
tend to keep the simple intuitive picture from leading order in mind, but one should be aware that 
the fundamental relationship is now more complicated, and that the parton densities no longer have a 
clear probabilistic interpretation. 

The reader should also be aware that there is currently a significant amount of debate concerning 
many aspects of PDF definitions and usage: 

- The "initial condition" for the PDFs, i.e., their shape in x at some low value of Q'j., and other 
constraints imposed on their evolution, such as positivity, flavour symmetries, treatment of mass 
effects, and extrapolation beyond the fit region. Each PDF group has its own particular ideol- 
ogy when it comes to these issues, and while the differences caused by these choices in well- 
constrained regions may appear small, the user should be warned that large differences can 
occur when extrapolating, e.g., to small x, or for observables that are particularly sensitive, e.g., 
to flavour symmetries, etc. 

- Using PDFs extracted using higher-order matrix elements in lower-order calculations, as, e.g., 
when using NLO PDFs as input to an LO calculation. In principle, the higher-order PDFs are 
better constrained and the difference between, e.g., an NLO and an LO set should formally be 
beyond LO precision, so that one might be tempted to simply use the highest-order available 
PDFs for any calculation. However, as described in section 2.4, it is often possible to partly ab- 
sorb higher-order terms into lower-order coefficients. In the context of PDFs, the fit parameters 
of lower-order PDFs will effectively attempt to "compensate" for missing higher-order contri- 
butions in the matrix elements. To the extent those higher-order contributions are universal, this 
is both desirable and self-consistent. However, this will only give an improvement when used 
with matrix elements at the same order as those used to extract the PDFs. It is therefore quite 
possible that NLO PDFs used in conjunction with LO matrix elements give a worse agreement 
with data than LO PDFs do. 




(48) 



Physical Observable 



Extracted Quantity 



25 



- PDF uncertainties. Uncertainty estimates for PDF determinations is a highly delicate procedure, 
owing in part to the diversity of the data sets that enter into the fitting procedures (especially 
since some data sets appear to have "tensions", i.e., mutual incompatibilities, between them), 
but also the differences in philosophy mentioned above (e.g., on parametrizations and evolution 
constraints) can cause apparent incompatibilities between different sets which are hard to give 
precise uncertainty estimates for. Currently, a consensus on meaningful uncertainty estimates is 
slowly building, though future years are likely to see continued active discussions on how best 
to address this topic. 

- How to use PDFs in conjunction with parton-shower Monte Carlo codes. The initial-state show- 
ers in a Monte Carlo model are essentially supposed to mimic the evolution in the PDFs, and 
vice versa. However, since PDF fits are not done with MC codes, but instead use analytical 
resummation models that are not identical to their MC counterparts, the PDF fits are essentially 
"tuned" to a slightly different resummation than that incorporated in a given MC model. Since 
both types of calculations are supposed to be accurate at least to LL, any difference between 
them should in principle be subleading. In practice, not much is known about the size and im- 
pact of this ambiguity, so we mention it mostly to make sure the reader is aware that it exists. 
Known differences include: the size of phase space (purely coUinear massless PDF evolution vs. 
the finite-transverse-momentum massive MC phase space), the treatment of momentum conser- 
vation and recoil effects, additional higher-order effects explicitly or implicitly included in the 
MC evolution, choice of renormalization scheme and scale, and, for those MC algorithms that 
do not rely on collinear (DGLAP, see [ ]) splitting kernels (e.g., the various kinds of dipole 
evolution algorithms, see [ ]), differences in the effective factorization scheme. 



3.2 Elastic and Inelastic Components of atot 
Elastic scattering consists of all reactions of the type 

A{pa)B{pb) ^ Aip'^)B{p's) , (49) 

where A and B are particles carrying momenta pA and ps, respectively. Specifically, the only ex- 
changed quantity is momentum; all quantum numbers and masses remain unaltered, and no new 
particles are produced. Inelastic scattering covers everything else, i.e., 

AB ^ X ^ AB , (50) 

where X ^ AB signifies that one or more quantum numbers are changed, and/or more particles are 
produced. The total hadron-hadron cross section can thus be written as a sum of these two physically 
distinguishable components, 

O-tot(s) = O-el(s) + Crinel(s) , (51) 

where s = {pA + Ps)^ is the beam-beam centre-of-mass energy squared. 

If A and/or B are not elementary, the inelastic final states may be further divided into "diffrac- 
tive" and "non-diffractive" topologies. This is a qualitative classification, usually based on whether 
the final state looks Uke the decay of an excitation of the beam particles (diffractive^"), or not (non- 
diffractive), or upon the presence of a large rapidity gap somewhere in the final state which would 
separate such excitations. 

'"An example of a process that would be labeled as diffractive would be if one the protons is excited to a which then 
decays back to +71"", without anything else happening in the event. In general, a whole tower of possible diffractive 
excitations are available, which in the continuum limit can be described by a mass spectrum falling roug hly as Ah'P /M^. 
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Given that an event has been labeled as diffractive, either within the context of a theoretical 
model, or by a final-state observable, we may distinguish between three different classes of diffractive 
topologies, which it is possible to distinguish between physically, at least in principle. In double- 
diffractive (DD) events, both of the beam particles are diffractively excited and hence none of them 
survive the collision intact. In single-diffractive (SD) events, only one of the beam particles gets 
excited and the other survives intact. The last diffractive topology is central diffraction (CD), in which 
both of the beam particles survive intact, leaving an excited system in the central region between them. 
(This latter topology includes "central exclusive production" where a single particle is produced in the 
central region.) That is, 

O-inel(s) = CrSD(s) + O-Dd(s) + O-Cd(s) + CJndIs) , (52) 

where "ND" (non-diffractive, here understood not to include elastic scattering) contains no gaps in 
the event consistent with the chosen definition of diffraction. Further, each of the diffractively ex- 
cited systems in the events labeled SD, DD, and CD, respectively, may in principle consist of several 
subsystems with gaps between them. Eq. (52) may thus be defined to be exact, within a specific 
definition of diffraction, even in the presence of multi-gap events. Note, however, that different theo- 
retical models almost always use different (model-dependent) definitions of diffraction, and therefore 
the individual components in one model are in general not directly comparable to those of another. 
It is therefore important that data be presented at the level of physical observables if unambiguous 
conclusions are to be drawn from them. 

3.3 Minimum-bias and soft inclusive physics 

The term "minimum-bias" (MB) is an experimental term, used to define a certain class of events that 
are selected with the minimum possible trigger bias, to ensure they are as inclusive as possible-'. In 
theoretical contexts, the term "minimum-bias" is often used with a slightly different meaning; to de- 
note specific (classes of) inclusive soft-QCD subprocesses in a given model. Since these two usages 
are not exactly identical, in these lectures we have chosen to reserve the term "minimum bias" to 
pertain strictly to definitions of experimental measurements, and instead use the term "soft inclusive" 
physics as a generic descriptor for the class of processes which generally dominate the various exper- 
imental "minimum-bias" measurements in theoretical models. This parallels the terminology used in 
the review [ ], from which most of the discussion here has been adapted. See equation (52) above for 
a compact overview of the types of physical processes that contribute to minimum-bias data samples. 
For a more detailed description of Monte Carlo models of this physics, in particular ones based on 
Multiple Parton Interactions (MPI), see Section 4.4. 

3.4 Underlying event and jet pedestals 

In events containing a hard parton-parton interaction, the underlying event (UE) can be roughly con- 
ceived of as the difference between QCD with and without including the remnants of the original 
beam hadrons. Without such "beam remnants", only the hard interaction itself, and its associated 

^'A typical min-bias trigger would thus be the requirement of at least one measured particle in a given rapidity region, 
so that all events which produce at least one observable particle would be included, which must, indeed, be considered 
the minimal possible bias. In principle, everything is a subset of minimum-bias, including both hard and soft processes. 
However, compared to the total minimum-bias cross section, the fraction that is made up of hard processes is only a very 
small tail. Since only a tiny fraction of the total minimum-bias rate can normally be stored, the minimum-bias sample would 
give quite poor statistics if used for hard physics studies. Instead, separate dedicated hard-process triggers are typically 
included in addition to the minimum-bias one, in order to ensure maximal statistics also for hard physics processes. 
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parton showers and hadronization, would contribute to the observed particle production. In reality, 
after the partons that participate in the hard interaction have been taken out, the remnants still contain 
whatever is left of the incoming beam hadrons, including also a partonic substructure, which leads to 
the possibility of "multiple parton interactions" (MPI), as will be discussed in section 4.4. Due to the 
simple fact that the remnants are not empty, an "underlying event" will always be there — but how 
much additional energy does it deposit in a given measurement region? A quantifation of this can be 
obtained, for instance, by comparing measurements of the UE to the average activity in minimum-bias 
events at the same t/s. Interestingly, it turns out that the underlying event is much more active, with 
larger fluctuations, than the average MB event. This is called the jet pedestal effect (hard jets sit on 
top of a higher-than-average "pedestal" of underlying activity), and is interpreted as follows. When 
two hadrons collide at non-zero impact parameter, high-px interactions can only take place inside 
the overlapping region. Imposing a hard trigger therefore statistically biases the event sample toward 
more central collisions, which will also have more underlying activity. See Section 4.4 for a more 
detailed description of Monte Carlo models of this physics, based on MPI. 

4 Monte Carlo Event Generators 

In this section, we discuss the physics of Monte Carlo generators and their mathematical foundations, 
at an introductory level. We shall attempt to convey the main ideas as clearly as possible without bury- 
ing them in an avalanche of technical details. References to more detailed discussions are included 
where applicable. We assume the reader is already familiar with the contents of the preceding sections 
of this report, in particular section 2.3 on matrix elements and section 2.4 on parton showers. Several 
of the discussions rely on material from the recent more comprehensive review in [ ] , which also 
contains brief descriptions of the physics implementations of each of the main general-purpose event 
generators on the market, together with a guide on how to use (and not use) generators in various 
connections, and a collection of comparisons to important experimental distributions. We highly rec- 
ommend readers to obtain a copy of that review, as it is the most comprehensive and up-to-date review 
of event generators currently available. Another useful and pedagogical review on event generators is 
contained in the 2006 ESHEP lectures by Sjostrand [42], with a more recent update in [45]. 

4.1 Perturbation Theory with Markov Chains 

Consider again the Born-level cross section for an arbitrary hard process, F, differentially in an arbi- 
trary infrared-safe observable O, as obtained from equation (34): 



where the integration runs over the full final-state on-shell phase space of F (this expression and 
those below would also apply to hadron collisions were we to include integrations over the parton 
distribution functions in the initial state), and the 6 function projects out a 1 -dimensional slice defined 
by O evaluated on the set of final-state momenta which we denote ^p. 

To make the connection to parton showers, we insert an operator, S, that acts on the Bom-level 
final state before the observable is evaluated, i.e.. 




Born 




(53) 




(54) 
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Formally, this operator — the evolution operator — will be responsible for generating all (real and 
virtual) higher-order corrections to the Bom-level expression. The measurement 5 function appearing 
explicitly in equation (53) is now implicit in S. 

Algorithmically, parton showers cast S as an iterative Markov (i.e., history-independent) chain, 
with an evolution parameter, Qe, that formally represents the factorization scale of the event, below 
which all structure is summed over inclusively. Depending on the particular choice of shower algo- 
rithm, Qe may be defined as a parton virtuality (virtuality-order showers), as a transverse-momentum 
scale (px-ordered showers), or as a combination of energies times angles (angular ordering). Regard- 
less of the specific form of Qe, the evolution parameter will go towards zero as the Markov chain 
develops, and the event structure will become more and more exclusively resolved. A transition from 
a perturbative evolution to a non-perturbative one can also be inserted, when the evolution reaches an 
appropriate scale, typically around 1 GeV. This scale thus represents the lowest perturbative scale that 
can appear in the calculations, with all perturbative corrections below it summed over inclusively. 

Working out the precise form that S must have in order to give the correct expansions discussed 
in section 2.4 takes a bit of algebra, and is beyond the scope we aim to cover in these lectures. Heuris- 
tically, the procedure is as follows. We noted that the singularity structure of QCD is universal and that 
at least its first few terms are known to us. We also saw that we could iterate that singularity structure, 
using universality and unitarity, thereby bootstrapping our way around the entire perturbative series. 
This was illustrated by the right-hand pane of figure 10 in section 2.4. 

Skipping intermediate steps, the form of the all-orders pure-shower Markov chain, for the evo- 
lution of an event between two scales Qei > Qe2, is, 

S{<!?f,Qei.Qe2,0) = A($p,(3bi,<3£2)*(0-0($p)) 



F + exclusive above Q e2 



E2 d1>^ 



Sri'^E+l) A($i., Qe1,Qf+i) S{<^f+1,Qf+1, Qe2, O) , 



F + 1 inclusive above Q e2 

(55) 



with the so-called Sudakov factor, 

A{<^p,Qei,Qe2) = exp 



r 



Qe2 



(56) 



defining the probability that there is no evolution (i.e., no emissions) between the scales Qei and 
Qe2, according to the radiation functions Sr to which we shall return below. The term on the first 
line of equation (55) thus represents all events that did not evolve as the resolution scale was lowered 
from Qei to Qe2, while the second line contains a sum and phase-space integral over those events 
that did evolve — including the insertion of 5(<I>p^_^) representing the possible further evolution of 
the event and completing the iterative definition of the Markov chain. 

The factor d^p_^^ /d^p defines the chosen phase space factorization. Our favourite is the so- 
called dipole- antenna factorization, whose principal virtue is that it is the simplest Lorentz invariant 
factorization which is simultaneously exact over all of phase space while only involving on-shell 
momenta. For completeness, its form is 

— , — = — ^ = dsai dsih 7, — , (57) 

d^p d^2 2^ WlT^Sr 
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Fig. 12: Illustration of the double-counting problem caused by naively adding cross sections involving 
matrix elements with different numbers of legs. 



which involves just one colour-anticolour pair for each r, with invariant mass squared Sr = {pa + 
Pi + Pb)'^- Other choices, such as purely collinear ones (only exact in the collinear limit or involving 
explicitly off-shell momenta), more global ones involving all partons in the event (more complicated, 
in our opinion), or less global ones with a single parton playing the dominant role as emitter, are also 
possible, again depending on the specific algorithm considered. 

The radiation functions Sr obviously play a crucial role in these equations, driving the emission 
probabilities. For example, if 5^ — )• 0, then A — exp(O) = 1 and all events stay in the top line of 
equation (55). Thus, in regions of phase space where Sr is small, there is little or no evolution. 
Conversely, for Sr — )• oo, we have A — > 0, implying that all events evolve. One possible choice for 
the radiation functions Sr was implicit in equation (40), in which we took them to include only the 
leading (double) singularities, with r representing colour-anticolour pairs. In general, the shower may 
exponentiate the entire set of universal singular terms, or only a subset of them (for example, the terms 
leading in the number of colours Nc), which is why we here let the explicit form of Sr be unspecified. 
Suffice it to say that in traditional parton showers, Sr would simply be the DGLAP splitting kernels 
(see, e.g., [22]), while they would be so-called dipole or antenna radiation functions in the various 
dipole -based approaches to QCD (see, e.g., [ , 35,46-49]). 

The procedure for how to technically "construct" a shower algorithm of this kind, using random 
numbers to generate scales distributed according to equation (55), is described more fully in [34], 
using a notation that closely parallels the one used here. The procedure is also described at a more 
technical level in the review [ ], though using a slightly different notation. Finally, a pedagogical 
introduction to Monte Carlo methods in general can be found in [_ v]. 

4.2 Matching 

The essential problem that leads to matrix-element/parton-shower matching can be illustrated in a 
very simple way. Assume we have computed the LO cross section for some process, F, and that we 
have added an LL shower to it, as in the left-hand pane of figure 12. We know that this only gives us an 
LL description of F + 1. We now wish to improve this from LL to LO by adding the actual LO matrix 
element for F + 1. Since we also want to be able to hadronize these events, etc, we again add an LL 
shower off them. However, since the matrix element for F + 1 is divergent, we must restrict it to cover 
only the phase-space region with at least one hard resolved jet, illustrated by the half-shaded boxes in 
the middle pane of figure 12. Adding these two samples, however, we end up counting the LL terms 
of the inclusive cross section for F + 1 twice, since we are now getting them once from the shower 
off F and once from the matrix element for F + 1, illustrated by the dark shaded (red) areas of the 
right-hand pane of figure 1 2. This double-counting problem would grow worse if we attempted to add 
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Fig. 13: Illustration of the original matctiing scheme implemented in HERWIG [51, 52], in which the 
dead zone of the Herwig shower was used as an effective "matching scale" for one emission beyond 
a basic hard process. 



more matrix elements, with more legs. The cause is very simple. Each such calculation corresponds 
to an inclusive cross section, and hence naive addition would give 

Ctot = '7o;incl + '7i;incl = C^O;excl + 2(Ti;incl • (58) 

Instead, we must match the coefficients calculated by the two parts of the full calculation — showers 
and matrix elements — more systematically, for each order in perturbation theory, so that the nesting 
of inclusive and exclusive cross sections is respected without overcounting. 

Given a parton shower and a matrix-element generator, there are fundamentally three different 
ways in which we can consider matching the two [ ] : 

1. Slicing: The most commonly encountered matching type is currently based on separating 
(slicing) phase space into two regions, one of which is supposed to be mainly described by hard matrix 
elements and the other of which is supposed to be described by the shower. This type of approach 
was first used in HERWIG [.">i ], to include matrix-element corrections for one emission beyond the 
basic hard process [^1,52]. This is illustrated in figure 13. The method has since been generalized 
by several independent groups to include arbitrary numbers of additional legs [53-57]. Effectively, 
the shower approximation is set to zero above some scale (either due to the presence of explicit dead 
zones in the shower, as in HERWIG, or by vetoing any emissions above a certain matching scale, as in 
the (L)-CKKW [53, 54, 56] and MLM [ , ] approaches), causing the matched result to be identical 
to the matrix element (ME) in that region, modulo higher-order corrections. We may sketch this as 

ME corrections 
Matched (above matching scale) = Exact x {l + 0{as)) , (59) 

where the "shower-corrections" include approximate Sudakov factors and reweighting factors ap- 
plied to the matrix elements in order to obtain a smooth transition to the shower-dominated region. 
Below the matching scale, the small difference between the matrix elements and the shower approx- 
imation can be dropped (since their leading singularities are identical and this region by definition 
includes no hard jets), yielding the pure shower answer in that region, 

shower correction 

Matched (below matching scale) = Approximate + (Exact — Approximate) 

= Approximate + non-singular 
— ^ Approximate . (60) 
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Fig. 14: Illustration of slicing approaches to matching, with up to two additional emissions beyond 
the basic process. The showers off F and F + 1 are set to zero above a specific "matching scale". 
(The number of coefficients shown was reduced a bit in these plots to make them fit in one row.) 
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Fig. 15: Illustration of the Mc@NLO approach to matching. In the middle pane, cyan boxes de- 
note non-singular correction terms, while the egg-coloured ones denote showers off such corrections, 
which cannot lead to double-counting at the LL level. 



This type of strategy is illustrated in figure 14. Since this strategy is discontinuous across phase 
space, a main point here is to ensure that the behaviour across the matching scale be as smooth as 
possible. CKKW showed [ ] that it is possible to remove any dependence on the matching scale 
through NLL precision by careful choices of all ingredients in the matching; technical details of the 
implementation (affecting the 0{as) terms in eq. (59)) are important, and the dependence on the 
unphysical matching scale may be larger than NLL unless the implementation matches the theoretical 
algorithm precisely [54, 56, 58]. One should also be aware that all strategies of this type are quite 
computing intensive. This is basically due to the fact that a separate phase-space generator is required 
for each of the n-parton correction terms, with each such sample a priori consisting of weighted events 
such that a separate unweighting step (often with quite low efficiency) is needed before an unweighted 
sample can be produced. 

2, Subtraction: Another way of matching two calculations is by subtracting one from the other 
and correcting by the difference, schematically 

shower correction 

' ^ ^ ^ 

Matched = Approximate + (Exact — Approximate) . (61) 

This looks very much like the structure of an NLO fixed-order calculation, in which the shower ap- 
proximation plays the role of subtraction terms, and indeed this is what is used in strategies like 
Mc@NLO [59-61], illustrated in figure 15. In this type of approach, however, negative-weight events 
will generally occur, for instance in phase-space points where the approximation is larger than the 
exact answer. This motivated the development of the so-called POWHEG approach [62], illustrated 
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Fig. 16: Illustration of the PowHEG approach to matching. In the middle pane, cyan boxes de- 
note non-singular correction terms, while the egg-coloured ones denote showers off such corrections, 
which cannot lead to double-counting at the LL level. 
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Fig. 17: Illustration of the MENLOPS approach to matching. Note that each of the POWHEG and 
CKKW samples are composed of separate sub-samples, as illustrated in figures 14 and 16. 



in figure 1 6, which is constructed specifically to prevent negative- weight events from occurring and 
simultaneously to be more independent of which parton-shower algorithm it is used with. The ad- 
vantage of these methods is obviously that NLO corrections to the Born level can be systematically 
incorporated. However, a systematic way of extending this strategy beyond the first additional emis- 
sion is not available, save for combining them with a slicing-based strategy for the additional legs, as 
in Menlops [63], illustrated in figure 17. These issues are, however, no more severe than in ordinary 
fixed-order NLO approaches, and hence they are not viewed as disadvantages if the point of reference 
is an NLO computation. 

3. Unitarity: The oldest, and in our view most attractive, approach [64,65] consists of working 
out the shower approximation to a given fixed order, and correcting the shower splitting functions at 
that order by a multiplicative factor given by the ratio of the matrix element to the shower approxima- 
tion, phase-space point by phase-space point. We may sketch this as 

correction 

shower , ^ ^ 

V Exact 

Matched = Approximate x . (62) 

Approximate 

When these correction factors are inserted back into the shower evolution, they guarantee that the 
shower evolution off n — 1 partons correctly reproduces the n-parton matrix elements, without the 
need to generate a separate n-parton sample. That is, the shower approximation is essentially used 
as a pre-weighted (stratified) all-orders phase-space generator, on which a more exact answer can 
subsequently be imprinted order by order in perturbation theory. In the original approach [64, 65], 
used by Pythia [66,67], this was only worked out for one additional emission beyond the basic hard 
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Fig. 18: Illustration of the two purely unitarity -based approaches to matching discussed in the text. 
Only one event sample is produced by each of these methods, and hence no sub-components are 
shown. 



process. In POWHEG [62], it was extended to include also virtual corrections to the Born-level matrix 
element. Finally, in ViNCiA [- '], it has been extended to include arbitrary numbers of emissions at 
tree level, though that method has so far only been applied to final-state showers. An illustration of the 
perturbative coefficients that can be included in each of these approaches is illustrated in figure 18, as 
usual with green (darker shaded) boxes representing exact coefficients and yellow (light shaded) boxes 
representing logarithmic approximations. Finally, two more properties unique to this method deserve 
mention. Firstly, since the corrections modify the actual shower evolution kernels, the corrections are 
automatically resummed in the Sudakov exponential, which should improve the logarithmic precision 
once A; > 2 is included, and secondly, since the shower is unitary, an initially unweighted sample of 
(n — l)-parton configurations remains unweighted, with no need for a separate event-unweighting or 
event-rejection step. 

4.3 The String Model of Hadronization 

In the context of event generators, hadronization denotes the process by which a set of post-shower 
partons is transformed into a set of primary hadrons, which may then subsequently decay further. This 
non-perturbative transition takes place at the hadronization scale, which by construction is identical 
to the infrared cutoff of the parton shower. In the absence of a first-principles solution to the relevant 
dynamics, event generators use QCD-inspired phenomenological models to describe this transition. 

Although non-perturbative QCD is not solved, we do have some knowledge of the properties 
that such a solution must have. For instance, Poincare invariance, unitarity, and causality are all 
concepts that apply beyond perturbation theory. In addition, lattice QCD provides us a means of 
making explicit quantitative studies in a genuinely non-perturbative setting, albeit only of certain 
questions. 

An important result in "quenched" lattice QCD^^ is that the potential of the colour dipole field 
between a charge and an anticharge appears to grow linearly with the separation of the charges, when 
the separation is greater than about a femtometer.This is known as "linear confinement", and it forms 
the starting point for the string model of hadronization. 

Starting from early concepts developed by Artru and Mennessier [68], several hadronization 
models based on strings were proposed in the late 1970'ies and early 80'ies. Of these, the most 
sophisticated and widely used today is the so-called Lund model, implemented in the Pythia code. 
We shall therefore concentrate on that particular model here, though many of the overall concepts 

^^Quenched QCD implies no "dynamical" quarks, i.e., no g ^ qq splittings allowed. 
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Fig. 19: Illustration of the transition between a Coulomb potential at short distances to the string-like 
one of equation (63) at large qq separations. 



would be shared by any string-inspired method. (A more extended discussion can be found in the 
very complete and pedagogical review of the Lund model by Andersson [69].) 

Consider the production of a qq pair from vacuum, for instance in the process e^e~ — 7*/Z — 
qq — )• hadrons. As the quarks move apart, hnear confinement implies that a potential 

V{r) = Kr (63) 

is asymptotically reached for large distances, r. (At short distances, there is a Coulomb term propor- 
tional to 1/r as well, but this is neglected in the Lund model.) This potential describes a string with 
tension (energy per unit length) k. The physical picture is that of a colour flux tube being stretched 
between the q and the q, figure 19. From hadron mass spectroscopy the string tension k, is known to 
be 

K ~ IGeV/fm ~ 0.2 GeV^. (64) 

A straightforward Lorentz-invariant description of this object is provided by the massless relativistic 
string in l-i-l dimensions, with no transverse degrees of freedom. The mathematical, one-dimensional 
string can be thought of as parameterizing the position of the axis of a cylindrically symmetric flux 
tube. (Note that the expression "massless" is somewhat of a misnomer, since k, effectively corresponds 
to a "mass density" along the string.) 

As the q and q move apart, their kinetic energy is gradually converted to potential energy, stored 
in the growing string spanned between them. In the "quenched" approximation, in which g ^ qq split- 
tings are not allowed, this process would continue until the endpoint quarks have lost all their momen- 
tum, at which point they would reverse direction and be accelerated by the now shrinking string. In 
the real world, quark-antiquark fluctuations inside the string field can make the transition to become 
real particles by absorbing energy from the string, thereby screening the original endpoint charges 
from each other and breaking the string into two separate colour-singlet pieces, (qq) — (qif) + (q'q), 
illustrated in figure 20 a. This process then continues until only ordinary hadrons remain. (We will 
give more details on the individual string breaks below.) More complicated multi-parton topologies 
including gluons are treated by representing gluons as transverse "kinks". Thus soft gluons effec- 
tively "build up" a transverse structure in the originally one-dimensional object, with infinitely soft 
ones absorbed into the string without leading to modifications. For strings with finite-energy kinks, 
the space-time evolution is then slightly more involved [(i*-)], and modifications to the fragmentation 
model to handle stepping across gluon corners have to be included, but the main point is that there 
are no separate free parameters for gluon jets. Differences with respect to quark fragmentation arise 
simply because quarks are only connected to a single string piece, while gluons have one on either 
side, increasing the energy loss per unit (invariant) time from a gluon to the string by a factor of 2 
relative to quarks, which can be compared to the ratio of colour Casimirs Ca/Cf = 2.25. 
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Fig. 20: a) Illustration of string breaking by quark pair creation in the string field, b) Illustration of 
the algorithmic choice to process the fragmentation from the outside-in, splitting off a single on-shell 
hadron in each step. 



Since the string breaks are causally disconnected (as can easily be realized from space-time 
diagrams [ ]), they do not have to be considered in any specific time-ordered sequence. In the 
Lund model, the string breaks are instead generated starting with the leading hadrons, containing the 
endpoint quarks, and iterating inwards towards the centre of the string, alternating randomly between 
fragmentation off the left- and right-hand sides, respectively, figure 20b. This has the advantage that a 
single on-shell hadron can be split off in each step, making it straightforward to ensure that only states 
consistent with the known spectrum of hadron resonances are produced, as will be discussed below. 

The details of the individual string breaks are not known from first principles. The Lund model 
invokes the idea of quantum mechanical tunneling, which leads to a Gaussian suppression of the 
energies and masses imparted to the produced quarks. 



Pioh{m pj_ ) oc exp exp , (65) 



where ruq is the mass of the produced quark and p± is the transverse momentum imparted to it by the 
breakup process (the antiquark obviously has the same mass and opposite p±). 

Due to the factorization of the p±a.nd m dependence implied by equation (65), the spectrum 
of produced quarks in this model is independent of the quark flavour, with a universal average value 
of 

(pig) = = k/tt {250MeYf . (66) 

Bear in mind that "transverse" is here defined with respect to the string axis. Thus, the p^in a frame 
where the string is moving is modified by a Lorentz boost factor. Also bear in mind that is here 
a purely non-perturbative parameter. In a Monte Carlo model with a fixed shower cutoff, the effec- 
tive amount of "non-perturbative" p±ma.y be larger than this, due to effects of additional unresolved 
soft-gluon radiation below the shower cutoff scale. In principle, the magnitude of this additional com- 
ponent should scale with the cutoff, but in practice it is up to the user to enforce this by retuning the 
effective a parameter when changing the hadronization scale. Since hadrons receive p± contributions 
from two breakups, one on either side, their average transverse momentum squared will be twice as 
large, 

(pi^> = 2a^ . (67) 
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The mass suppression implied by equation (65) is less straightforward to interpret. Since quark 
masses are notoriously difficult to define for light quarks, the value of the strangeness suppression 
must effectively be extracted from experimental measurements, e.g., of the K/tt ratio, with a resulting 
suppression of roughly s/u ~ s/d ~ 0.2 - 0.3. Inserting even comparatively low values for the 
charm quark mass in equation (65), however, one obtains a relative suppression of charm of the order 
of 10"^^. Heavy quarks can therefore safely be considered to be produced only in the perturbative 
stages and not by the soft fragmentation. 

Baryon production can be incorporated in the same basic picture [~ by allowing string breaks 
to occur also by the production of pairs of so-called diquarks, loosely bound states of two quarks in 
an overall 3 representation (e.g., red + blue = antigreen). Again, the relative rate of diquark-to-quark 
production is not known a priori and must be extracted from experimental measurements, e.g., of the 
p/vr ratio. More advanced scenarios for baryon production have also been proposed, in particular the 
so-called popcorn model [7 1 , 72], which is normally used in addition to the diquark picture and then 
acts to decrease the correlations among neighbouring baryon-antibaryon pairs by allowing mesons to 
be formed inbetween them. Within the Pythia framework, a fragmentation model including explicit 
string junctions [ ] has so far only been applied to baryon-number- violating new -physics processes 
and to the description of beam remnants (and then acts to increase baryon stopping [74]). 

This brings us to the next step of the algorithm, assignment of the produced quarks within 
hadron multiplets. The fragmenting quark (antiquark) may combine with the antiquark (quark) from a 
newly created breakup to produce either a vector or a pseudoscalar meson, or, if diquarks are involved, 
either a spin- 1/2 or spin-3/2 baryon. Unfortunately, the string model is entirely unpredictive in this 
respect, and this is therefore the sector that contains the largest amount of free parameters. From 
spin counting alone, one would expect the ratio V/S of vectors to pseudoscalars to be 3, but in 
practice this is only approximately true for B* /B. For lighter flavours, the difference in phase space 
caused by the V-S mass splittings implies a suppression of vector production. Thus, for D* / D, the 
effective ratio is already reduced to about ~ 1.0 - 2.0, while for K* /K and p/tt, extracted values 
range from 0.3 - 1.0. Recall, as always, that these are production ratios of primary hadrons, hence 
feed-down complicates the extraction of these parameters from experimental data, in particular for 
the lighter hadron species. The production of higher meson resonances is assumed to be low in a 
string framework^\ For diquarks, separate parameters control the relative rates of spin-1 diquarks vs. 
spin-0 ones and, likewise, have to extracted from data, with resulting values of order {qq)i/{qq)o ~ 
0.075-0.15. 

With p'j_ and now fixed, the final step is to select the fraction, z, of the fragmenting end- 
point quark's longitudinal momentum that is carried by the created hadron. In this respect, the string 
picture is substantially more predictive than for the flavour selection. Firstly, the requirement that the 
fragmentation be independent of the sequence in which breakups are considered (causality) imposes 
a "left-right symmetry" on the possible form of the fragmentation function, f{z), with the solution 



which is known as the Lund symmetric fragmentation function (normalized to unit integral). As a 
by-product, the probability distribution in invariant time r of q'q breakup vertices, or equivalently 
r = (kt)'^, is also obtained, with dP/dV oc exp(— ftF) implying an area law for the colour flux, 
and the average breakup time lying along a hyperbola of constant invariant time tq ~ lO^^'^s [69]. 

^^The four L = 1 multiplets are implemented in PYTHIA, but are disabled by default, largely because several states are 
poorly known and thus may result in a worse overall description when included. 




(68) 
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Fig. 21 : Illustration of the iterative selection of flavours and momenta in the Lund string fragmentation 
model. 



The a and b parameters are the only free parameters of the fragmentation function, though a may 
in principle be flavour-dependent. Note that the explicit mass dependence in f{z) implies a harder 
fragmentation function for heavier hadrons (in the rest frame of the string). 

The iterative selection of flavours, p± , and z values is illustrated in figure 2 1 . A parton produced 
in a hard process at some high scale Q\jy emerges from the parton shower, at the hadronization scale 
QiR, with 3-momentum p = {p±Q,p^), where the "+" on the third component denotes "light-cone" 
momentum, p± = E ± p^. Next, an adjacent dd pair from the vacuum is created, with relative 
transverse momenta ±p±i. The fragmenting quark combines with the d from the breakup to form a 
TT"*", which carries off a fraction zi of the total lightcone momentum The next hadron carries off a 
fraction Z2 of the remaining momentum, etc. 

For massive endpoints (e.g., c and b quarks, or hypothetical hadronizing new -physics particles, 
generally called i?-hadrons), which do not move along straight lightcone sections, the exponential 
suppression with string area leads to modifications of the form [75], f{z) — )• f{z)/z^"^Q, with mg 
the mass of the heavy quark. Strictly speaking, this is the only fragmentation function that is consistent 
with causality in the string model, though a few alternative forms are typically provided as well. 

Note, however, that the t&na fragmentation function in the context of non-perturbative hadroniza- 
tion models is used to denote only the corrections originating from scales below the infrared cutoff 
scale of the parton shower. That is, the fragmentation functions introduced here are defined at an 
intrinsically low scale of order Q ~ 1 GeV. It would therefore be highly inconsistent and misleading 
to compare them directly to those that are used in fixed-order and/or analytical-resummation contexts, 
which are typically defined at a factorization scale of order the scale of the hard process. 

4.4 Multiple Parton Interactions 

In Monte Carlo modeling contexts, multiple parton interactions (MPI) denote the possibility of hav- 
ing multiple partonic 2 — )• 2 interactions occurring within a single hadron-hadron collision. The 
most striking and easily identifiable consequence of MPI is thus arguably the possibility of observing 
several distinct (i.e., hard) parton-parton interactions in some fraction of hadron-hadron events. Ad- 
ditional jet pairs produced in this way are sometimes referred to as "minijets", but in the interest of 
maintaining a compact terminology, we shall here just call them MPI jets. The main distinguishing 
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feature of such jets is that they tend to form back-to-back pairs, with httle total p±. For comparison, 
jets from bremsstrahlung tend to be aligned with the direction of their "parent" partons. The fraction 
of multiple interactions that give rise to additional reconstructible jets is, however, quite small (how 
small depends on the exact jet definition used). Additional soft interactions, below the jet cutoff, are 
much more plentiful, and can give significant corrections to the colour flow and total scattered en- 
ergy of the event. This affects the final-state activity in a more global way, increasing multiplicity 
and summed Et distributions, and contributing to the break-up of the beam remnant in the forward 
direction. 

The first detailed Monte Carlo model for perturbative MPI was proposed by Sjostrand and van 
Zijl in [76], and with some variation this still forms the basis for most modern implementations. Here, 
we therefore focus on that model and on its more recent "interleaved" version [ ] . Some discussion 
of alternative models as well as additional references to the history and development of the subject of 
multiple interactions can be found in [ . ]. 

An intuitive way of arriving at the idea of multiple interactions is to view hadrons simply as 
'bunches' of incoming partons. No physical law then prevents several distinct pairs of partons from 
undergoing scattering processes within one and the same hadron-hadron collision. The other key idea 
to bear in mind is that the exchanged QCD particles are coloured, and hence such multiple interactions 
— even when soft — can cause non-trivial changes to the colour topology of the colliding system as 
a whole, with potentially major consequences for the particle multiplicity in the final state. 

To begin to construct a model for this, we first observe that, at low p±, t-channel propagators 
almost go on shell (reminiscent of the case of bremsstrahlung, described in detail in section 2.4), 
which causes the differential QCD parton-parton scattering cross sections (such as the Rutherford one 
illustrated in section 2) to become very large, behaving roughly as: 

di dp^. 

An integration of this cross section from a lower cutoff p_Lmin to ^/s, using the full (leading- 
order) QCD 2 — 2 matrix elements folded with some recent parton-density sets, is shown in figure 22, 
for pp collisions at 14 TeV [ ]. The solid curves, representing the calculated cross sections as 
functions of p_Lmin> are compared to a few different predictions for fitot (the total pp cross section 
[79]), shown as horizontal lines with different dashing styles on the same plot. Physically, the jet 
cross section can of course not exceed the total pp one, yet this is what appears to be happening at 
scales of order 4-5 GeV in figure 22. How to interpret this behaviour? 

Recall that the interaction cross section is an inclusive number. Thus, if a single hadron-hadron 
event contains two parton-parton interactions, it will count twice in (J2^2 but only once in cJtot> and 
so on for higher parton-parton interaction numbers. In the limit that all the individual parton-parton 
interactions are independent and equivalent (to be improved on below), we have 

0-2^2(P±min) = {n-){p±min) <7tot , (70) 

with (n) (p_Lmin) giving the average of a Poisson distribution in the number of parton-parton interac- 
tions above pxmin per hadron-hadron collision, 

-n / \ 1/ \f Mn exp [-(n)(p_Lmin)] 

Vn[P±min) = [{n){p±ram)] , , (71) 

n] 

and that number may well be above unity. This simple argument in fact expresses unitarity; instead 
of the total interaction cross section diverging as p_Lmin (which would violate unitarity), we have 
restated the problem so that it is now the number of interactions per collision that diverges. 
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Fig. 22: The inclusive jet cross section calculated at LO for three different proton PDFs, compared 
to various extrapolations of the non-perturbative fits to the total pp cross section at 14 TeV centre-of- 
mass energy. From [78]. 



Two important ingredients remain to be introduced in order to fully regulate the remaining 
divergence. Firstly, the interactions cannot use up more momentum than is available in the parent 
hadron. This will suppress the large-n tail of the naive estimate above. Obviously, exact momentum 
conservation is included in all Monte Carlo models currently on the market, although the details vary 
somewhat from model to model. In the PYTHlA-based models [^6, 77, 8' ], the multiple interactions 
are ordered mp±, and the parton distributions for each successive interaction are explicitly constructed 
so that the sum of x fractions can never be greater than unity. In the Herwig models [81,82], instead 
the uncorrected estimate of (n) above is used directly as an initial guess, but the actual generation of 
interactions stop once the energy-momentum conservation limit is exceeded (with the last "offending" 
interaction also removed from consideration). 

The second ingredient suppressing the number of interactions, at low p± and x, is colour screen- 
ing; if the wavelength ~ of an exchanged coloured parton becomes larger than a typical colour- 
anticolour separation distance, it will only see an average colour charge that vanishes in the limit 
p± — )• 0, hence leading to suppressed interactions. This screening effectively provides an infrared 
cutoff for MPI similar to that provided by the hadronization scale for parton showers. A first estimate 
of an effective lower cutoff due to colour screening would be the proton size 

h 0.2GeV-fm , 

P±min - - ~ — TT^r ~ 0.3 GeV ~ Aqcd , (72) 

Tp U.7 im 

but empirically this appears to be too low. In current models, one replaces the proton radius rp in 
the above formula by a "typical colour screening distance" d, i.e. an average size of a region within 
which the net compensation of a given colour charge occurs. This number is not known from first 
principles, so effectively this is simply a cutoff parameter, which can then just as well be put in 
transverse momentum space. The simplest choice is to introduce a step function Q{p± —p±mm)^ such 
that the perturbative cross section completely vanishes below the p±mm scale. Alternatively, one may 
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note that the jet cross section is divergent like a'^{p'j)/p'^, cf. eq. (69), and that therefore a factor 

would smoothly regularize the divergences, now with p±Q as the free parameter to be tuned to data. 
Regardless of whether it is imposed as a smooth (Pythia and Sherpa) or steep (Herwig ++) func- 
tion, this is one of the main "tuning" parameters in such models. Note also that this parameter does not 
have to be energy-independent. Higher energies imply that parton densities can be probed at smaller 
X values, where the number of partons rapidly increases. Partons then become closer packed and the 
colour screening distance d decreases. The uncertainty on the energy and/or x scaling of the cutoff is 
a major concern when extrapolating between different collider energies. 

We now turn to the origin of the so-called "pedestal effect", the observational fact that hard jets 
appear to sit on top of a higher "pedestal" of underlying activity than events with no hard jets. This 
is interpreted as a consequence of impact-parameter-dependence, as follows. In peripheral collisions, 
only a small fraction of events contain any high-p^ activity, whereas central collisions are more likely 
to contain at least one hard scattering; a high-px triggered sample will therefore be biased towards 
small impact parameters. The ability of a model to describe the shape of the pedestal (e.g., to describe 
both "minimum-bias" data and underlying-event distributions simultaneously) therefore depends upon 
its modeling of the impact-parameter dependence, and correspondingly the impact-parameter shape 
constitutes another main tuning parameter for models that include this dependence. 

For each impact parameter, b, the number of interactions n can then still be assumed to be dis- 
tributed according to a Poissonian, eq. (71), again modulo momentum conservation, but now with the 
mean value of the Poisson distribution depending on impact parameter, If the matter distri- 

bution has a tail to infinity (as, e.g., Gaussians do), one may nominally obtain events with arbitrarily 
large b values. In order to obtain finite total cross sections, it is therefore necessary to give a separate 
interpretation to the "zero bin" of the Poisson distribution, which corresponds to no-interaction events. 
In models that attempt to describe the entire inelastic non-diffractive cross section, this bin is simply 
ignored, since the events in it can only represent elastic or diffractive scatterings, which are modeled 
separately. Alternatively, in models that pertain only to hard inelastic events, it can be reinterpreted as 
containing that fraction of the total inelastic cross section which do not contain any hard interactions. 

Finally, we should mention that there are two perturbative modeling aspects which go beyond 
the introduction of MPI themselves. In particular, this concerns 

1 . Parton showers off the MPI. 

2. Perturbative parton-rescattering effects. 

Without showers, MPI models would generate very sharp peaks for back-to-back MPI jets, 
caused by unshowered partons passed directly to the hadronization model. However, with the excep- 
tion of the oldest Pythia 6 model [ ], all of the general-purpose event- generator models do include 
such showers, and hence should exhibit more realistic (i.e., broader and more decorrelated) MPI jets. 
On the initial state side of the MPI shower issue, the main questions are whether and how correlated 
multi-parton densities are taken into account, and, as discussed previously, how the showers are regu- 
lated at low p± and/or low x. Although none of the MC models currently impose a rigorous correlated 
multi-parton evolution, all of them include some elementary aspects. The most significant for parton- 
level results is arguably momentum conservation, which is enforced explicitly in all the models. The 
so-called "interleaved" models [74,77] attempt to go a step further, generating an explicidy correlated 
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multi-parton evolution in which flavour sum rules can be imposed to conserve, e.g., the total numbers 
of valence and sea quarks across interaction chains. 

Perturbative rescattering in the final state occurs if partons are allowed to undergo several 
distinct interactions, with showering activity possibly taking place inbetween. This has so far not 
been studied extensively, but a first fairly complete model and exploratory study has been presented 
in the context of Pythia 8 [80]. In the initial state, parton rescattering effects have so far not been 
included in any of the general-purpose Monte Carlo models. 

4.5 Colour (Re)-Connections and Beam Remnants 

Consider now a hadron-hadron collision, i.e., including MPI, at the parton level, equivalent to a res- 
olution scale of about one GeV. The system of coloured partons emerging from the short-distance 
phase (primary parton-parton interaction plus parton-level underlying event plus beam-remnant par- 
tons) must now undergo the transition to colourless hadrons. Infrared sensitive observables, such as 
individual hadron multiplicities and spectra are crucially dependent on the parton-parton correlations 
in colour space, and on the properties and parameters of the hadronization model used. Here, we 
concentrate on the specific issues connected with the structure of the event in colour space. 

Keeping the short-distance parts unchanged, the colour structure inside each of the MPI systems 
is normally still described using just the ordinary leading-colour matrix-element and parton-shower 
machinery described in sections 2.3 and 2.4. The crucial question, in the context of MPI, is then how 
colour is neutrahzed between different MPI systems, including also the remnants. Since these systems 
can lie at very different rapidities (the extreme case being the two opposite beam remnants), the strings 
spanned between them can have very large invariant masses (though normally low p_[_), and give rise 
to large amounts of (soft) particle production. Indeed, in the context of soft-inclusive physics, it is 
precisely these "inter-system" strings which furnish the dominant particle production mechanism, and 
hence their modeling is an essential part of the infrared physics description. 

As discussed more fully in [41], there is a large amount of ambiguity concerning how to ad- 
dress this, and a substantial amount of variation between current models. Experimental investigations 
of colour reconnections at LEP [ - ] were only able to exclude some fairly extreme models, with 
comparatively moderate ones still allowed. Furthermore, in hadron collisions the initial state con- 
tains soft colour fields with wavelengths of order the confinement scale. The presence of such fields, 
unconstrained by LEP measurements, could impact in a non-trivial way the process of colour neu- 
tralization [87, 88]. And finally, the MPI produce an additional amount of displaced colour charges, 
translating to a larger density of hadronizing systems. It is not known to what extent the collective 
hadronization of such a system differs from a simple sum of independent systems. 

A new generation of colour-reconnection toy models have therefore been developed specifically 
with soft-inclusive and underlying-event physics in mind [ - ], and also the cluster-based [ ] and 
Generalized- Area-Law [93] models have been revisited in that context. Although still quite crude, 
these models do appear to be able to describe significant features of the Tevatron and LHC data, such 
as the (p±)(A'ch) distribution in minimum-bias data, which appears to be quite sensitive to this effect. 
It is nonetheless clear that the details of the full fragmentation process in hadron-hadron collisions are 
still far from completely understood. 

4.6 Tuning 

The main virtue of general-purpose Monte Carlo event generators is their ability to provide a complete 
and fully differential picture of collider final states, down to the level of individual particles. This 
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allows them to be used as detailed — albeit approximate — theoretical references for measurements 
performed at accelerators like the LHC, against which models of both known and 'new' physics can 
be tested. As has been emphasized in these lectures, the achievable accuracy depends both on the 
inclusiveness of the chosen observable and on the sophistication of the simulation itself. An important 
driver for the latter is obviously the development of improved theoretical models, e.g., by including 
matching to higher-order matrix elements, more accurate resummations, or better non-perturbative 
models, as discussed in the previous sections; but it also depends crucially on the available constraints 
on the remaining free parameters of the model. Using existing data to constrain these is referred to as 
generator tuning. 

Although Monte Carlo models may appear to have a bewildering array of independently ad- 
justable parameters, it is worth keeping at the front of one's mind that most of these parameters only 
control relatively small (exclusive) details of the event generation. The majority of the (inclusive) 
physics is determined by only a few, very important ones, such as, e.g., the value of the strong cou- 
pling, in the perturbative domain, and the form of the fragmentation function for massless partons, in 
the non-perturbative one. 

Armed with a good understanding of the underlying model, an expert would therefore normally 
take a highly factorized approach to constraining the parameters, first constraining the perturbative 
ones and thereafter the non-perturbative ones, each ordered in a measure of their relative significance 
to the overall modeling. This factorization, and carefully chosen experimental distributions corre- 
sponding to each step, allows one to concentrate on just a few parameters and distributions at a time, 
reducing the full parameter space to manageable-sized chunks. Still, each step will often involve more 
than one single parameter, and non-factorizable corrections still imply that changes made in subse- 
quent steps can change the agreement obtained in previous ones by a non-negligible amount, requiring 
additional iterations from the beginning to properly tune the entire generator framework. 

Recent years have seen the emergence of automated tools that attempt to reduce the amount of 
both computer and manpower required for this task, for instance by making full generator runs only 
for a limited set of parameter points, and then interpolating between these to obtain approximations 
to what the true generator result would have been for any intermediate parameter point, as, e.g., in 
the Professor tool [ , ]. Automating the human expert input is of course more difficult. In the 
tools currently on the market, this question is addressed by a combination of input solicited from 
the generator authors (e.g., which parameters and ranges to consider, which observables constitute a 
complete set, etc) and the elaborate construction of non-trivial weighting functions that determine how 
much weight is assigned to each individual bin and to each distribution. The field is still burgeoning, 
however, and future sophistications ai^e to be expected. Nevertheless, at this point the overall quality 
of the tunes obtained with automated methods appear to at least be competitive with the manual ones. 

A sketch of a reasonably complete tuning procedure, without going into details about the pa- 
rameters that control each of these sectors in individual Monte Carlo models, would be the following: 

1) Keep in mind that inabilities of models to describe data is a vital part of the feedback cycle 
between theory and experiment. Also keep in mind that perturbation theory at LOxLL is doing very 
well if it gets within 10% of a given IR safe measurement. An agreement of 5% should be considered 
the absolute sanity limit, beyond which it does not make any sense whatsoever to tune further. The 
advent of NLO Monte Carlos may reduce these numbers slightly, but only for quantities for which 
one expects NLO precision to hold. However, the sanity limit should be taken to be at least twice as 
large for quantities governed by non-perturbative physics. For some quantities, e.g., ones for which 
the underlying modeling is known to be poor, an order-of-magnitude agreement or worse may have 
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to be accepted. Attempting to force Monte Carlo models to describe data far outside their domains 
of validity must be expected to produce similar side effects as attempting to turn a Fiat into a Ferrari 
merely by cranking up the engine revolutions. 

2) Final-State Radiation and Hadronization: mainly using LEP and other e^e~ collider data. 
On the IR safe side, there are event shapes and jet observables, the latter including rates, resolutions, 
masses, shapes, and jet-jet correlations. On the IR sensitive side, special attention should be paid to 
the high-z tail of the fragmentation spectra, where a single hadron carries a large fraction of an entire 
jet's momentum, since this is the tail that is most likely to give "fake jets". Depending on the focus 
of the tuning, attention should also be paid to identified-particle rates and ratios, and to fragmentation 
in events containing heavy quarks and/or gluon jets. Usually, more weight is given to those particles 
that are most copiously produced, though this again depends on the focus. Finally, particle-particle 
correlations and baryon production are typically some of the least well constrained components of the 
overall modeling. The scaling properties of IR safe vs. IR sensitive contributions can be tested by 
comparing data at several different e+e~ collider energies. 

3) Initial-State Radiation, and so-called "Primordial""* kT": here, one would in principle 
like to use data from DIS reactions, which are less complicated to interpret than full hadron-hadron 
collisions. However, due to difficulties in translating between the ep and pp environments, this is 
normally not what is done in practice. Instead, the main constraining distribution is the dilepton 

distribution in Drell-Yan events in hadron-hadron collisions. For any observables containing ex- 
plicit jets, be aware that the underlying event can produce small horizontal shifts in jet p± distribu- 
tions, which may in turn result in seemingly larger-than-expected vertical changes if the distributions 
are falling sharply. Also note that the ISR evolution is sensitive to the choice of PDFs, with caveats as 
discussed in section 3.1. 

4) Initial-Final Connections: e.g., radiation from colour lines connected to the initial state 
and jet broadening in hadron collider environments. This is one of the most poorly controlled parts 
of most MC models. Keep in mind that it is not directly constrained by pure final-state observables, 
such as LEP fragmentation, nor by pure initial-state observables, such as the Drell-Yan p± spectrum, 
which is why we list it as a separate item here. In principle, DIS would again be a prime territory for 
placing constraints on this aspect at least for quark jets, but in practice more often inclusive-jet and 
other multi-jet processes (such as jets) in hadron colliders are used. 

5) Underlying Event: Good constraints on the overall level of the underlying event can be 
obtained by counting the summed transverse energy (more IR safe) and/or particle multiplicities and 
average transverse momenta (more IR sensitive) in regions transverse to a hard trigger jet (more IR 
safe) or particle (more IR sensitive). Constraints on the fluctuations of the underlying event are also 
important, and can be obtained, e.g., by comparing to measurements of the RMS of such distributions. 
Again, note that the UE is sensitive to the choice of PDFs. 

6) Colour (Re-)Connections and other Final-State Interactions: By Final-State Interactions, 
we intend a broad spectrum of possible collective effects that may be included to a greater or lesser 
extent in various models. These effects include: Bose-Einstein correlations, colour reconnections, 
hydrodynamics, string interactions, Cronin effect, etc. As a rule, these effects are non-perturbative 
and hence should not modify IR safe observables appreciably. They can, however, have drastic effects 
on IR sensitive ones, such as particle multiplicities, and particle momentum distributions, wherefore 
useful constraints are typically furnished by particle-particle correlations, by measurements of particle 

Primordial kr'- an additional soft p± component that is injected on top of the p± generated by the initial-state shower 
itself, see [-n, Section 7.1]. 
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momentum spectra as functions of quantities believed to serve as indicators of the strength of these 
phenomena (such as event multiplicity), and/or by collective-flow-type measurements. Finally, if the 
model includes a universal description of underlying event and soft-inclusive QCD, as many MPI- 
based models do, then minimum-bias data can also be used as a control sample, though one must 
then be careful either to address diffractive contributions properly or to include only data samples that 
minimize their impact. 

7) Beam Remnants: Constraints on beam remnant fragmentation are most easily obtained in 
the forward region, but, e.g., the amount of baryon transport from the remnant to a given rapidity re- 
gion can also be used to probe how much the colour structure of the remnant was effectively disturbed, 
with more baryon transport indicating a larger amount of "beam baryon blowup". 

We round off by emphasizing that comparisons of specific models and tunes to data can be 
useful both as immediate tests of commonly used models, and to illustrate the current amount of 
theoretical uncertainty surrounding a particular distr ibution. Independently of how well the models fit 
the data, such comparisons also provide a set of well-defined theoretical reference curves that serve 
as useful guidelines for future studies. However, the conclusions that can be drawn from comparisons 
of individual tunes of specific models on single distributions are necessarily limited. In order to 
obtain more general conclusions, a strategy for a more coherent and over-arching look at both the 
data and the models was recently proposed in [ ]. Specifically, rather than performing one global 
tune to all the data, as is usually done, a more systematic check on the validity of the underlying 
physics model can be obtained by instead performing several independent optimizations of the model 
parameters for a range of different phase space windows and/or collider environments. In regions in 
which consistent parameter sets are obtained, with predictions that are acceptably close to the data, the 
underlying model can then be considered as interpolating well, i.e., it is universal. If not, a breakdown 
in the ability of the model ability to span different physical regimes has been identified, and can be 
addressed, with the nature of the deviations giving clues as to the nature of the breakdown. With 
the advent of automated tools making it easier to run several optimizations without much additional 
computing overhead, such systematic studies are now becoming feasible, with a first example given 
in [96]. 
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